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1.  Objectives. 

The  main  goal  of  the  Project  is  developing  a  solid  theoretical  and  numerical  background  for  the 
analysis  of  control  systems  with  spatially  distributed  feedback.  The  subjects  of  particular  interest  are 
the  study  of  new  approaches  to  classical  Fourier  filtering  and  irreversible  spatial  transforms  as  novel 
types  of  distributed  feedback  control.  For  such  a  type  of  analysis,  new  effective  numerical  methods  for 
processing  information  in  the  frequency  domain  should  be  developed.  These  numerical  techniques  can 
be  useful  in  multi-frame  alignment  for  low-light  applications  (optical  coherence  tomography),  stereo¬ 
image  analysis  for  robotics,  automated  medical  image  processing,  texture  matching,  etc.  Therefore,  the 
objectives  split  into  three  closely  connected  tasks,  accomplishment  of  which  could  provide  substantial 
improvement  of  the  performance  and  robustness  of  distributed  feedback-controlling  systems  subjected 
to  micro-control  devices. 

Task  1:  Distributed  Fourier-Hermite  Projection  Filtering. 

Nowadays,  “texture”  is  an  important  tenn  in  the  world  of  computer  vision,  graphical 
applications,  computer  modeling,  and  optical  system  controlling.  The  task  of  texture  processing  first 
came  up  with  the  task  of  making  automatic  image  recognition  systems.  Textures  in  the  context  of 
image  recognition  can  be  used  in  different  ways:  image  segmentation;  detecting  the  physical  properties 
and  materials  of  the  surfaces  represented  in  the  image;  detecting  the  shape  and  position  of  the  objects 
and  surfaces  in  the  image.  These  problems  need  also  an  effective  method  of  edge  detection  that  is 
stable  to  noise. 

Our  main  goal  was  to  develop  a  new  effective  method  for  the  tasks  of  image  filtering,  texture 
identification  and  edge  detection. 

The  first  approach  to  implement  the  above  goals  was  based  on  Hermite  series  expansion  time- 
frequency  filtering  techniques.  Non-hierarchical  and  hierarchical  texture  coding  schemes  using  this 
approach  have  been  developed.  This  approach  allowed  us  also  to  create  a  new  stable  method  of  edge 
detection. 

The  second  approach  is  based  on  Tikhonov’s  regularization  method.  New  image  filtering 
algorithms  using  analytically  found  Green’s  functions  have  been  created  and  implemented.  A 
numerical  method  for  solving  the  corresponding  Euler-Lagrange  equation  has  been  also  designed. 

Task  2:  Controlling  Discrete  Fourier  Filtering. 

The  development  of  a  new  class  of  optical  filters  based  on  Fourier  series  with  a  discrete  filtering 
kernel  (discrete  filters)  instead  of  traditional  integral  convolution  attracts  a  substantial  practical 
interest.  Discrete  spatial  filters  are  much  easier  to  fabricate  using  photolithography  technology. 
Contemporary  technology  allows  producing  filters  with  a  certain  number  of  “steps”  for  a  micron  of 
height.  On  the  other  hand,  there  are  technological  capabilities  for  producing  high-speed  dynamic 
discrete  Fourier  filters,  using  MEMS  devices.  Consequently,  the  possibility  to  use  filters  with  a 
discrete  number  of  steps  in  thickness  and  a  finite  number  of  elements  in  the  transverse  plane  is  a 
favorable  option,  and  such  filters  should  be  modeled  numerically.  In  the  framework  of  our  Project, 
spatial  filters  are  expected  to  enhance  the  characteristics  of  2D  nonlinear  optical  systems  with  a 


1 


controlling  feedback  loop.  Thus,  specialized  software  is  required  to  model  the  dynamics  of  such  a 
system  containing  a  spatial  filter  in  its  feedback  loop  as  well  as  the  theoretical  background  of  the 
methods  implemented. 

Task  3:  Controlling  Irreversible  Transforms  of  Spatial  Arguments. 

Large-scale  spatial  transformations  in  the  feedback  loop  create  novel  possibilities  for  synthesizing 
light  fields  with  desired  properties,  as  compared  with  conventional  spatial  filtering  techniques. 
Contraction,  reflection,  shift,  and  rotation  are  examples  of  reversible  large-scale  transforms  capable  of 
dramatically  changing  the  spatio-temporal  dynamics  of  the  feedback  system  and  actually  appear  as  a 
novel  and  convenient  controlling  tool.  More  advantages  could  be  accomplished  by  applying 
irreversible  transforms  resulting  in  the  formation  of  localized  and  periodic  ID  optical  patterns.  A  new 
mathematical  approach  to  give  an  adequate  description  of  irreversible  transforms  was  recently 
suggested  by  participants  of  the  Project.  The  approach  is  based  on  the  theory  of  distributions  and  takes 
into  account  a  wide  range  of  Lebesgue-measurable  2D  irreversible  transforms  including  focusing  onto 
a  point  or  a  planar  curve.  It  is  thus  necessary  to  accomplish  comprehensive  investigation  of  optimal 
control  problems  for  optical  systems  governed  by  PDEs  with  controlling  reversible  and  irreversible  2- 
D  argument  transforms  in  the  feedback  loop.  Consequently,  a  solid  mathematical  background  and 
specialized  numerical  methods  are  required  to  study  the  corresponding  optimization  algorithms  and  to 
apply  them  to  optical  information  processing. 

2.  Technical  Description  of  Work  Accomplished 
2.1.  Distributed  Fourier-Hermite  Projection  Filtering 

In  the  scope  of  the  Project,  non-hierarchical  and  hierarchical  texture  coding  schemes  using 
Hermite  series  expansion  time-frequency  approach  have  been  developed.  This  approach  allowed  us  to 
create  a  new  stable  method  for  edge  detection.  A  new  method  of  image  filtering  and  edge  detection 
based  on  the  Tikhonov  regularization  method  has  been  also  designed. 

1.0.0.  Non-hierarchical  texture  coding  scheme 


A  filtering  method  for  texture  analysis  has  been  elaborated.  It  is  based  on  applying  a  set  of 
band-pass  filters  to  extract  the  so-called  feature  vectors  from  the  image.  Each  filter  responds  most 
strongly  to  a  selected  spatial-frequency  and  orientation  bands.  In  our  work,  Hermite  expansion  series 
are  taken  to  define  the  filters.  The  Hermite  functions  are  defined  as: 

y/n  (x)  =  ^ ^  e  -  ■  —  - — - .  They  can  be  also  calculated  recursively: 

V2  "nl.Jn  dx 


Vo  = 


■e 


-x1!  1 


Vx  = 


-x  12 


Vn  =  XJ~  ■  Vn-X  -  J—  ■  Vn-2  ,  ^  2 

V  n  V  n 

1-D  Hermite  functions,  y/n  (jc,  y)  =  y/n  (x)  •  1 ,  have  a  strongly  pronounced  orientation.  To  perform  an 

adequate  analysis,  the  functions  should  be  rotated  by  certain  angles.  During  the  first  stage  of 
investigations,  it  was  found  that  8  orientation  angles,  0,  22.5,  45,  67.5,  90,  112.5,  135,  and  157.5 
degrees,  are  enough  to  obtain  stable  representation  of  texture  features. 

In  this  approach,  to  get  the  feature  vectors,  we  considered  functions  i//„(x,y),  where  n= 0..64,  and 
introduced  6  energy  coefficients: 

Ei  =  (ao)"  +  (ai)2 , 

E2  =  (<X2)“+(a3)“  +  ((X4)2, 


2 


E3-  (as)-  +(a6)2+(ay)“+(a8)2, 

Ee  =  (a33)“+(<X34)~  +  •  •  ■  +  (a63)“  +(ct64)2  , 


where  ( x,y )  is  the  source  image,  a,  are  the  Hermite  series  coefficients.  Hence,  with  the  8  orientations 
we  get  a  48-dimensional  feature  vector  for  each  texture  example.  An  analysis  shows  (see  Fig.  1.1)  that 
the  feature  vectors  obtained  by  this  filtering  algorithm  are  suitable  for  texture  parameterization. 


Fig.  1.1  Two  examples  of  the  same  texture  and  corresponding  feature  vectors.  In  the  2-D 
representations  of  the  feature  vectors,  spatial  frequency  levels  increase  bottom-up,  and  the  orientation 

changes  from  left  to  right,  being  wrapped  around. 


2.0.0.  Hierarchical  texture  coding  scheme 

When  decomposing  a  function  by  the  standard  scheme  and  using  the  coefficients  of  1-D 
Hermite  series  expansion  without  rescaling,  some  problems  arise  with  texture  analysis.  The  most 
significant  of  these  is  that  images  taken  from  the  real  world  are  unstructured,  i.e.  the  luminosity,  scale, 
size,  and/or  perspective  of  objects  may  vary,  causing  inaccuracy  during  computer  processing.  To  avoid 
the  above-mentioned  problem,  we  used  the  idea  of  hierarchical  coding. 

This  approach  allows  us  to  use  several  sets  of  patterns,  which  correspond  to  different  features 
the  image  may  contain.  These  sets  differ  in  spatial  scale.  1-D  Hermite  series  expansion  with  a 
changeable  spatial  scale  enables  solving  these  problems. 

To  compare  texture  parameterization  methods,  a  special  program  module  has  been  developed. 
It  allows  loading  grayscale  texture  samples  and  applying  to  them  three  different  texture 
parameterization  methods.  As  a  result,  the  program  returns  a  table  of  comparison  of  the  methods 
applied  to  the  loaded  texture  patterns.  An  example  of  such  a  comparison  table  is  shown  in  Fig.  1.2 


a 

B 

Ordinary  Coding 

Hierarchical  Coding 

Hierarchical  Coding  without  subtractions 

At 

A2 

0.004505 

0.005349 

0.003739 

B1 

B2 

0.009905 

0.008160 

0.007742 

Cl 

C2 

0.017778 

0.011848 

0.009946 

A3 

A4 

0.008807 

0.006047 

0.002422 

B3 

B4 

0.020075 

0.010667 

0.006275 

C3 

C4 

0.128322 

0.101591 

0.080176 
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A1 

B1 

0.488420 

0.492238 

0.491653 

B1 

Cl 

0.315644 

0.304597 

0.305442 

Fig.  1.2  Differences  in  processing  some  texture  samples  with  different  coding  schemes.  Brodatz 


textures  of  the  same  type  are  denoted  by  the  same  letters. 


The  comparison  table  contains  the  absolute  values  of  the  differences  between  the  feature 
vectors  of  the  loaded  texture  patterns.  The  difference  is  calculated  using  the  following  formula: 

48 

£(<*,- A)2 

p(a,b)  = - - — - ,  where  a,  and  /?,  are  the  Hermite  series  coefficients  for  the 

max 

1=1  i=i 

corresponding  coding  scheme. 

The  analysis  of  the  comparison  table  shows  that  the  best  method  to  be  applied  to  texture 
parameterization  and  segmentation  tasks  is  the  method  called  “Hierarchical  Coding  without 
subtractions”.  Assuming  this,  we  have  chosen  this  method  to  solve  the  texture  segmentation  problem. 

To  perfonn  the  texture  segmentation  task  using  the  “Hierarchical  Coding  without  subtractions” 
method,  another  program  module  has  been  developed.  It  allows  one  to  load  a  grayscale  picture, 
fracture  it  into  blocks  of  a  selected  size,  tune  the  method  parameters,  and  perform  the  segmentation. 
Some  results  of  the  segmentation  program  operation  are  shown  in  Fig.  1.3  and  Fig.  1.4 


Fig.  1.3.  The  results  of  applying  the  texture  segmentation 
algorithm  with  different  block  texture  differences  thresholds 
to  a  Brodatz  texture  patterns-based  image. 
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Fig.  1.4.  The  results  of  applying  the  texture  segmentation 
algorithm  with  different  blocks  texture  differences  thresholds  to  a 
natural-scene  image. 


Hermite  function-based  parameterization  implies  some  calculating  problems.  The  main 
problem  is  that  the  calculation  of  the  Hermite  function  Hn(x)  value  at  a  given  point,  x,  is  very  machine¬ 
time  expensive  for  a  large  n.  So,  we  were  forced  to  calculate  these  recursively.  To  solve  the  problem 
and  to  increase  the  calculation  performance  of  the  program,  we  developed  a  caching  calculation 
scheme,  which  is  based  on  the  fact  that  many  values  of  Hn(x)  are  retrieved  repeatedly  at  different 
stages  of  calculation.  Such  an  approach  requires  more  memory,  but,  on  the  other  hand,  it  allows  us  to 
increase  the  calculation  speed  up  to  10  times. 
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3.0.0.  Stable  edge  detection  methods 

The  problem  of  texture  analysis,  parameterization,  and  discrimination  using  1-D  and  2-D  Hermite 
series,  which  we  are  developing,  needs  a  stable  method  to  perform  image  edge  detection.  In  (J.F. 
Canny.  A  computational  approach  to  edge  detection  //  IEEE  Trans,  on  Pattern  Analysis  and  Machine 
Intelligence,  8(6):679— 698,  November  1986)  three  criteria  for  edge  detection  are  given.  These  have 
been  used  in  a  similar  form  in  a  good  deal  of  vision  research.  They  are:  1.  Good  detection.  There 
should  be  the  minimum  number  of  false  negatives  and  false  positives.  2.  Good  localization.  The  edge 
location  must  be  reported  as  close  as  possible  to  the  correct  position.  3.  Only  one  response  to  a  single 
edge.  Basically,  these  criteria  are  connected  with  the  mathematical  problem  of  stable  second  derivative 
zero  crossing  detection  for  a  perturbed  function.  Our  research  was  focused  on  developing  new 
algorithms  for  solving  this  problem  on  the  basis  of  the  Tikhonov  regularization  method  and  1-D 
Hermite  series  approach. 

A.  Edge  detection  method  based  on  Tikhonov’s  regularization 


The  method  based  on  Tikhonov’s  regularization  looks  as  follows:  for  a  function  y  e  Hl ,  we 
have  an  approximation  || y  -  ys ||  2  <  5  .  We  want  to  find  y  e  Hl ,  such  that  ||  v  -  y\H i  — »  0  for  S  -a  0. 
This  problem  is  being  solved  by  minimizing  the  Tikhonov  functional: 

>  l|2 


I  y-ysh 


+  a 


dx 


y 


— >  mm . 


The  problem  of  the  Tikhonov  functional  minimization  is  reduced  to  solving  the  Euler  equation  for  y 
(below  we  use  notation  ‘y’  instead  ): 

,  „  i  i 

\y  — y  =  — ys 

a  a 

v'(l)  =  y'(-l)  =  0 


|  -ay  +y  =  ys 
l/(l)  =  /(-l)  =  0 


Its  solution  is  ya(x)  =  J Ga(x,s)ys(s)ds,  with  the  kernel 


-1 


Ga(x,s)  = 


-fash 


— -ch— -  l)ch— + 1)  s<x 

2  Jfc  Jfc 


- fash 


-ch~=(s  -  l)c/?— + 1)  s>x 


The  graphs  of  the  kernel  look  as  follows: 
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1  ,0 


0,5 


a  -  0.1 


a-  10 


Fig.  1.5  Graphs  of  the  kernel  Ga(x,s),  - 1  <  jc  <  1,  - 1  <  5  <  1 . 

To  find  a  =  a(d) ,  the  solution  ya  (x)  must  satisfy  the  following  equation: 

1 

cp{a)  =  J (ya(x)  -  ya{x)fdx  =  S2  =>  a(5) 

-1 

The  proposed  method  finds  function  y  with  the  minimum  norm  of  the  derivative  in  L2  among  all  the 
functions  y(x)  from  the  set  ||y  -  y§  ||  2  ^  S . 


Fig.  1.6  Parameter  a  selection  by  discrepancy  method  ( S2  =0. 165,  a  =0.000986). 


A  very  interesting  consequence  of  the  above  approach  is  the  fact  that  the  solution  satisfies  the  equation 
-  ay"  +  y=ys-  So,  zero  crossing  of  the  second  derivative  is  equivalent  to  the  change  of  the  sign  of 

y  -  ys .  This  allows  us  to  construct  a  method  for  stable  second  derivative  zero  crossing  detection.  Its 

application  to  image  processing  gives  us  a  new  method  of  stable  edge  detection.  For  edge  detection  we 
use  Laplacian  zero-crossing  criterion.  From  the  equation 

-au"  +  u  =  us 


it  follows  that 


u"  =  0  u  -  us  =  0 . 
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So,  we  smooth  every  row  of  the  noisy  image  us ,  and  find  the  second  derivative  in  the  y  direction 

dru 


— y  .  Then,  we  smooth  every  column  of  the  initial  image  us  ,  and  find  the  second  derivative  in  the  x 


direction, 


d~u 

a 


.  Summing  these  results,  we  obtain  a  non-directional  second  derivative  (the  Laplacian 


operator:  A u  =  ).  Figure  1 .7  shows  applying  this  algorithm  to  edge  detection. 

dx  dy 


Fig.  1.7  Edge  detection  results. 


The  derivatives  used  for  edge  thresholding  can  be  calculated  by  the  following  formula: 
u  '(x)  = - - - [  sh  — (x  - 1  )ch  — (s  +  1)  Vg  ( s)ds - - - [  ch  — -!=  (s  - 1  )sh  ~^=  (x  + 1)  vs  ( s)ds . 


4a  sh 


4ash 


4a 


4a 


The  proposed  method  can  be  used  for  the  considered  task  of  edge  detection,  but  the  design  of 
2D  methods  can  make  the  algorithms  more  effective.  To  implement  this  approach,  variational 
regularization  methods  for  2D  image  filtering  were  considered.  Typical  variational  methods  of  image 
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restoration  obtain  a  filtered  version  of  some  degraded  image  us  as  the  minimizer,  ua  ,  of: 

Ea(u )  =  J [{u  -  ug )2  +  avF(|Vw|~ )]dx  . 
n 

The  first  summand  encourages  the  similarity  between  the  restored  image  ua  and  the  original  one  us, 

while  the  second  summand  rewards  smoothness.  The  smoothness  weight  a  >  0  is  called  the 
regularization  parameter.  In  our  case,  the  regularizer  T  is  supposed  to  satisfy  the  following  conditions: 

1.  T(.)  is  continuous  for  any  compact  K  c;  [0,oo). 

2.  XP(  | .  |2)  :  91"  — >•  5?  —  is  a  convex  function  . 

3.  T(.)  is  increasing  in  [0,oo). 

4.  There  exists  a  constant  ^>0:VF(52)>^52. 

For  this  class  of  regularization  methods,  one  can  establish  a  well-posed  and  scale-space  framework,  as 
for  nonlinear  diffusion  filtering.  The  regularization  parameter  a  should  be  considered  as  a  scale. 

Theorem  1,1,  (Properties  of  regularization  methods) 

(a)  (Well-posedness  and  regularity)  Let  us  e  L  00  (Q) .  Then  the  functional  Ea(u)  has  a  unique 

minimizer  ua  in  the  Sobolev  space  H 1  (Q) .  Moreover,  ua  e//2( Q)  and  |wa|i2(n)  continuously 
depends  on  a. 

(b)  (Extremum  principle)  Let  a  =  infa  us  and  b  =  supD  us.  Then  a<ua  <  b  on  Q. 

(c)  (Average  grey  level  invariance)  The  average  grey  level  ju  =  —  j '  us(x)dx  remains  constant  under 

Ml  n 

1  r 

regularization:  —  J  ua  (x)dx  =  /u. 

Ml  n 

(d)  (Convergence  to  a  constant  image  for  a  00  )  lima^.00|Ma  -  //||£P(n)  =  0  for  any  p  e  [l,oo). 

If  T  is  differentiable,  then  the  minimizer  of  Ea(u )  satisfies  the  Euler-Lagrange  equation: 


— 'A-E  =  div(W(\Vu\2)'Vu). 
a 

Summary  of  some  ¥  functions: 


Name  of  the  function  ‘P 

'P(t) 

T'(0 

Convexity 

Tikhonov 

T 

1 

Yes 

Total  Variation 

■ft 

1 

Yes 

2  ft 

Bouman  and  Sauer 

ta  0.5<a<l 

at"'1  0.5<a<l 

Yes 

Perona  and  Malik 

l-exp(-t) 

exp(-t) 

No 

Geman  and  McClure 

t 2 

1 

No 

l  +  t2 

(i  +  O2 
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Hebert  and  Leahy 

Log(l+t) 

2 

1  +  t 

No 

Green 

Log[ch(  yfi )] 

th(4t) 

Yes 

2  4t 

Hyper  surfaces 

Vi  +  f-i 

1 

Yes 

2Vl  +  t 

The  above  methods  were  numerically  analyzed.  The  Tikhonov  and  Perona-Malik  approaches  were 
found  as  the  most  promising.  The  first  of  these  approaches  (edge  detection  by  Tikhonov’s 
regularization)  is  a  2D  generalization  of  our  previous  ID  method.  In  this  case,  if  we  use  T'(t)=t,  the 
Euler-Lagrange  equation  is  reduced  to:  u-us  =  c/Au.  Therefore,  Aw  =  0  <=>  u  -us  =0.  The  proposed 
method  allows  us  to  obtain  edges  in  an  image  on  the  basis  of  Laplacian  zero  crossing  analysis. 


B.  Edge  detection  based  on  1-D  Hermite  series  approach 

The  Hermite  expansion  of  function /(A)  looks  as  follows: 

GO  00 

f(x)  =  Y,ai  •  Yffr)  ,  ai=  \  Wi(x)  ■  f(x)dx  , 

i= 0  —co 


where  a;  is  the  i- th  Fourier  coefficient,  and  %(x)  is  the  z'-th  Hermite  function.  Subsequently, 

f  OO  \  00 


f'(x) 


\k= 0  J  k= 0 


Basing  on  the  following  formula  that  is  valid  for  the  Hermite  functions 

¥n"  +  {2-n  +  l-x2)-¥„  =  0, 

00 

we  can  evaluate  /"(x)  —  ’  (-Y  ~  2  •  k  -1)  •  (x) . 

k-0 

Let  us  now  have  a  two-dimensional  image  I(x,y)  of  size  m  *  n.  Basing  on  the  above  formula  for 
one-dimensional  functions  decompositions,  we  define  the  following  scheme  for  processing  the  image 
I(x,  y).  The  rows  {L:  (x)}i=|  n  and  columns  {C;  m  are  apparently  analyzed  for  second  derivate 

zeros.  For  each  row  F,  (x)  we  have  a  set  of  second  derivative  values,  L”(x),  x  being  a  discrete 
argument:  x=l..m.  Having  the  set  of  second  derivative  values,  we  define  the  following  criterion  for  a 
zero-crossing  of  L'\x)  :  x  is  considered  a  zero-crossing  point  if  it  matches  the  inequality 

L .  ”  (x)  •  L"  (x  - 1)  <  -A  A  being  a  threshold  parameter. 

Second  derivative  zero  points  for  {Z.(x)}.=1  n  and  {C7  ( v)}  m  can  be  considered  as  edge  points  of  the 

source  image  I{x,  y).  Correspondingly  these  points  are  then  plotted  in  the  resulting  image  that  depicts 
the  edges  found  in  the  source  image. 

We  have  compared  the  described  method  of  Hennite  decomposition-based  edge  detection  with 
the  method  based  on  Fourier  decomposition.  The  method  differs  from  the  described  one  only  in 
obtaining  the  second  derivative  values: 
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/(*) 


^  (ak  cos  Ax  +  bk  sin  kx) , 

k=i 


2/o 

VTtt 


—  f  f(x)dx ,  ak  =  =  —  f  /(x)cos  Axt/x ,  bk  -  —  —  [ /(x)sin  kxdx  ■ 

7T  Jjr  Jl  _  Jrr  71 


Then  the  second  derivative  can  be  calculated: 


/"(*)  = 


—  +  cos  kx  +  bfr  sin  Ax)  "=  '^k2(-a/i  cos  kx-  \  sin  Ax) . 

,  2  A=1  J  A=1 


The  algorithm  of  finding  edges  using  this  formula  is  the  same  as  that  described  above  for  Hermite 
decomposition-based  edge  detection. 

Even  before  practical  comparison  of  using  the  Hennite  and  Fourier  series  for  edge  detection  in 
noisy  images,  we  can  see  a  great  advantage  of  Hermite  series  approach.  It  lies  in  the  fact  that  the 
coefficient  ak  of  the  function  /  expansion  is  multiplied  by  k  for  the  Hermite  series,  while  the  k 2 


multiplier  arises  instead  in  the  Fourier  series  approach.  This  means  that  the  Hermite  series  approach  is 
more  stable  when  dealing  with  noisy  data. 

While  having  almost  identical  results  when  detecting  edges  in  noise-free  images,  our  method 
displays  preferable  results  when  detecting  edges  in  noisy  images.  The  followings  Figures  illustrate  the 
results  obtained. 
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Fig.  1.8  Four  images  are  plotted:  source  image  of  a  white  circle,  with  a  20%  Gaussian  noise;  the 
results  of  190-harmonic  Fourier-based  edge  detection  method  application  (two  samples  with  different 
thresholds);  and  the  result  of  190  Hermite  functions  based  edge  detection  method  application.  It  can  be 
seen  that  the  Hermite  functions  based  method  gives  an  adequate  result  while  the  Fourier  method  does 
not  yield  adequate  edge  detection.  It  is  also  illustrated  that  changing  the  threshold  does  not  help  the 
Fourier  method  to  produce  good  results,  because  either  (with  a  small  threshold)  we  obtain  a  too  noisy 
resultant  image  that  does  not  emphasize  edges  or  (with  a  bigger  threshold)  we  obtain  an  image  with 
very  little  points  that  depict  the  edges  but  still  many  points  that  are  just  noise. 


Fig.  1.9  Source  image  “Lena”;  the  result  of  190-harmonic  Fourier-based  edge  detection  method 
application;  and  the  result  of  190  Hermite  functions  based  edge  detection  method  application. 


4.0.0.  Variational  methods  of  image  filtering 


The  proposed  method  based  on  Tikhonov’s  regularization  was  also  aimed  at  solving  the  edge 
detection  problem.  A  more  general  but  still  very  important  problem  of  image  filtering  can  be  solved  by 
minimizing  the  Tikhonov  functional  of  a  higher  order: 


Ea(u)=  Ti-us 


+  J3 


d- 


dx~ 


The  first  term  of  the  functional  is  a  data  fidelity  term,  while  the  second  rewards  smoothness  of  the 
solution.  This  problem  can  be  reduced  to  solving  the  Euler  equation  for  u  : 

4 


d4) 


+  4  a4u  =  u< 


w'(l)  =  =  0  ,  with  the  kernel  Ga(x,s)  = 

w'"(l)  =  w'"(-l)  =  0 


^a((5)M;(x)  x<s 

i= 1 
4 

YJbi(s)ui(x)  x>s 


i= 1 


The  following  procedure  has  been  developed  to  find  the  analytical  expression  for  the  kernel.  The 
fundamental  solutions  look  as  follows: 


ux  (x)  =  ea( x+1)  sin  a(x  +  1),  u2  (x)  =  ea{ x+1)  cos  a(x  +  1), 

w3(x)  =  e_a(v_1>  sina(x  - 1),  w4(x)  =  e~a(x~l)  cos a(x  -  1) 

After  substituting  4  =  ce(x  +  1)  rj  =  a(x  - 1) ,  their  derivatives  are  given  by 


u[(x)  =  ae*[  sin  4  + cos q] 

u"(x)  =  2 a2e^[  cosq] 

u'”(x)  =  2a 3  A  [-sin  4  +  cos 4] 

u'2(x)  =  ae ■"  [-  sin  4  +  cos  4] 

u"(x)  =  2«2e^[-sinq] 

u"Xx)  =  2«3e^[-sin4-cos4] 

u\  (x)  =  ae^  [-  sin  +  cos  r/\ 

u"(x)  =  2a1  A7  [-cos  77] 

m"(x)  =  2aie~v[  sin 7  + cos 7]] 

u\(x)  =  ae  ^f-sin^  -  cos 7] 

u\{x)  =  2a1e-\  sin  77] 

wj"(x)  =  2a  V “''[-sin^  +  cos^]] 
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Here,  c,  (5)  =  bj  (5) 


(s)  .  From  the  conjunction  conditions 


4 

^c,.(s)Mi(s)  =  °, 

1=1 

YJci(s)u’’(s)  =  0, 

.  1=1 


Y.CiisXis)  =  0, 


(=1 

4 


YjCjisy^s)  =  1 


i=l 


we  get 


g-a(s+l) 

Cj(s)  = - —  [-sina(s  + 1)  +  cosa(s  + 1)] 

8  a 

g-a(s+X) 

c2  (s)  = - — [-sina(s  +  l)-cosa(s  +  l)] 

8a 

c3  (5)  = - —  [sin  a(s  - 1)  +  cos  a{s  - 1)] 

8a 

ea(s- 1) 

c4(s )  = - —  [-sina(5-l)  +  cosa(5-l)] 

8a ' 

After  substituting  Ga(x,s)  in  the  boundary  conditions 
(G'a(-l,s)  =  G'a(l,S)  =  0, 

=  Gma(l,s)  =  0, 

we  finally  obtain  the  coefficients  bi(s),ai(s)  : 

p  =  c{(s)  +  e2u[c3(s)cos2a  +  c4(s)sin2a],  q  =  c2(s)  +  e2a[c3{s)sm2a  -  c4(s)  cos2a] 

,  ,  .  e4cTneos4a  +  nsin4al  -  n  ,  ,  .  e4cTnsin4a  -  q  cos  4a]  +  q 

b,(s)  = - t - - - r - — ,  bi(s)  = - - - 7 - - - 7 - 

l-2e4“  cos4a  +  e8“  “  1  -  2e4“  cos  4a  +  e8“ 

b3{s)  =  ela[-bx(s)cos2a  +  b2(s)sin2a],  b4(s)  =  e2a\  ^(sjsk^a +  b2(s)cos2a] 

We  are  flexible  to  change  the  boundary  conditions  for  the  Euler  equation.  For  the  case  of  the  following 
boundary  conditions 

fGB(-U)  =  GB(U)  =  0 

{G:(-l,s)  =  G:(l,s)  =  0, 

the  coefficients  h,  (5),  at  (5)  are  given  by  expressions: 

p  =  cx(s)- e2a[c3(s)cos2a  +  c4(s)sm2a],  q  =  c2(s)  -  e2u [c3(s)sm2a  -  c4(.s')cos2a] 

,  ,  .  e4cTneos4a  +  «sin4al  -  p  ,  ,  .  e4"  in  sin  4a  -  q  cos  4a]  +  q 

b,(s)  = - 7 - - - 7 - — ,  b2(s)  = - - - 7 - - - 7 - 

1  -2e4a  cos  4a  +  e8a  “  1  -  2e4“  cos  4a  +  e%a 

b3(s)  =  ela[bx(s)cos2a  -h2 (5) sin 2a],  b4(s )  =  -ela[bl{s)s’m2a +  b2{s)cos2a] 

The  above  versions  of  the  kernel  are  illustrated  in  Fig.  1.10. 
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Ga(x,s) ,  a  -  1 ,  boundary  conditions  with  even 
derivatives 


Ga  (x,  s) ,  a-  1 ,  boundary  conditions  with  odd 
derivatives 


0.10 

008 
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0.04 

0.02 

0.00 


X  Axis 


Fig.  1.10  Green’s  functions 


The  obtained  analytical  expressions  of  Green’s  functions  give  us  a  new  tool  for  solving  the  problem  of 
image  filtering. 

The  general  scheme  of  the  methods  based  on  Tikhonov’s  regularization  using  both  first  and 
second  order  terms,  which  we  analyzed  separately  before,  can  be  posed  as  a  problem  of  minimizing  the 
Tikhonov  functional: 


Ea(u)  =  f  ~us\\L2  +  « 


d_ 

dx 


+  J3 


L2 


d~ 


dx~ 


L2 


The  method  is  intended  to  be  also  used  in  the  edge  detection  procedure  based  on  the  analysis  of 
the  second  order  zero  crossings  of  the  intensity  function.  It  can  be  expected  that  the  method  based  on 
minimizing  the  first  derivative  ( /?  =0  )  will  suppress  noise  in  the  image  and  the  number  of  the  second 
order  zero  crossings  will  be  reduced.  The  method  with  the  second  derivative  ( a  =0  )  is  expected  to 
perform  stronger  smoothing,  while  preserving  the  second  order  zero  crossings. 


Here  we  can  consider  not  only  self-adjoint  statements  as  it  was  done  for  analytical 
investigations.  A  special  program  unit  has  been  designed  to  solve  this  problem.  We  will  describe  the 
algorithm  for  the  case 

I  J3i/4)  -  cai"  +  u  =  u s , 

\ u'{\)  =  u’(- 1)  =  0,  u"{  1)  =  w"(-l)  =  0  ' 


The  numerical  implementation  looks  as  follows.  By  approximating  the  derivatives  contained  in  the 
equation  with  finite  differences  and  proceeding  to  the  vector  notation  u,=u(ih )  (where  h  is  a  step  size, 
z'=0. . ./?),  we  can  rewrite  the  equation  as 

B  oc 

\(n  .  -?jy  4-  11  i  —  9  (u  -?jv  4- u  )+  (w.  -  2uM  +  W;+2  )]  +  — 


-2km  +k,)-2(km  -2 Ui+UM 
h 


i+2 


)]  +  TT  k-1  “2  U,+  Ui+ 1  ]  +  Ui=fi. 

h 


The  boundary  conditions  can  be  rewritten  as 

u0  -  U\  =0,  u0  -  2u1  +u2  =0,  un  -  un_x  =  0,  un  -  2zzn_j  +  un_2  =  0  . 

These  boundary  conditions  can  be  reduced  to  u0  =ux  =  u2  un  =  un_x  =  un  2 .  By  substituting  these 
conditions  into  the  equation,  we  obtain  an  n-2 -order  system  of  linear  equations 
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4-[3«2  -4m3  +u4]  +  7t\tu2  +u,]+u2  =f2 

n  n 

P 


—  [-3 u2  +6 «3  -4u4  +  i/5 ] h — -[u2  -  2u3  +  w4]  +  w3  =  /3 
/z4  /z 


13  cc 

7T k-2  -  4«W  +  6m,  -  4m!+1  +  M,.+2  ]  +  -y  [mm  -  2iii  +  uM  ]  +  ui=fi  i  =  4, 4 
n  n 


A[-3un-2  +6m„_3  -4 w„„4  +m„_5]  +  ^[m„_2  -2«„_3  +m„_4]  +  «„_3  = /„ 
h  h 

^[3m„_2  -4m„_3  +m„_4]  +  -^[-w„_2  +m„_3]  +  m„_2  =/„_2 


-3 


This  system  can  be  rewritten  in  a  compact  matrix  form  as  Au  -  f  . 

Further  analysis  has  been  performed  to  compare  the  influence  of  the  first  and  the  second  order 
term  in  the  Tikhonov  functional.  The  scheme  of  the  analysis  was  as  follows.  For  a  set  of  test  images, 
we  added  a  unifonnly  distributed  noise  (an  example  of  noisy  data  is  given  below): 


Fig.  1.11  Noisy  test  image. 

Then,  we  analyzed  the  value  of  the  PSNR  metric  between  the  initial  image  and  the  smoothed  image 
obtained  as  a  function  of  the  regularization  parameter  a  (for  the  case  of  the  pure  first  order  Tikhonov 
functional)  or  the  regularization  parameter  (3  (for  the  case  of  the  pure  second  order  Tikhonov 
functional).  It  was  found  that  the  obtained  maximal  values  of  PSNR  are  very  close  in  all  the  cases. 


Fig.  1.12  PSNR  values  as  functions  of  regularization  parameters. 


14 


Visual  analysis  of  the  best-PSNR  results  (Fig.  1.13)  was  perfonned  at  the  final  stage  to  compare 
features  of  the  first  and  second  order  Tikhonov  regularization  schemes. 


Fig.  1.13  First  order  (left)  and  second  order  (right)  Tikhonov  regularization  filtering  results. 

The  analysis  showed  that  the  first  order  Tikhonov  regularization  filtering  works  better  for  lower 
frequencies,  while  the  second  order  Tikhonov  regularization  filtering  is  preferable  to  treat  high- 
frequency  details. 
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1.0.  Controlling  Discrete  Fourier  Filtering. 

2.2.1.  Model  of  discrete  Fourier  filtering  in  nonlinear  optical  systems 

The  theoretical  background  of  optimal  Fourier  filtering  has  been  developed  and  applied  to  the 
problem  of  producing  phase-amplitude  Fourier  filters  aimed  at  forming  a  desired  phase  distribution  in 
a  light  wave  by  means  of  a  nonlinear  optical  system  with  such  a  Fourier  filter  in  the  feedback  loop 
(see.  Fig. 2.0). 


Fig.  2.0.  Conceptual  scheme  of  a  nonlinear  optical  system  with  a  spatial  filter  in  the  feedback  loop. 

Here,  a  spatial  filter,  SF,  is  placed  in  the  Fourier  plane  of  a  «4 -f»  system  that  produces  a  feedback 
signal  to  drive  a  nonlinear  medium,  NL.  In  the  framework  of  a  single -pass  approximation,  the 
dynamics  of  the  above  system  is  governed  by  a  Debye-type  nonlinear  parabolic  PDE.  This  equation 
governs  the  additional  phase  modulation,  u  =  u(x,  t),  of  a  light  wave  passed  through  a  thin  layer  of  a 
Kerr-type  nonlinear  medium  being  a  part  of  an  optical  feedback  system  with  a  Fourier  filter  in  its 
feedback  loop: 


t  —  +  u-IjAju  =  ®(u,T),  te[0,fi],  u(x,t  =  0)  =  w0(x), 
dt 


®  (u,T)  =  KSl 


A in 


+2  K8> 


2  Re(  A, At  {An  exp(z'w)})  +  KS3  Ft  {Ain  exp(z'w)} 


(2.1) 


Here,  r  is  the  characteristic  response  time  of  the  nonlinear  medium;  ld  is  the  diffusion  length  of  the 
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medium;  =  Vj_  is  the  transverse  Laplace  operator  with  respect  to  transverse  spatial  variable 
x  =  (xi,x2)  e  Q ,  where  Q  is  the  optical  system  aperture,  being  a  bounded  convex  2D  domain. 

Equation  (2.1)  is  supplemented  with  homogeneous  boundary  conditions  of  the  Dirichlet, 
Neumann,  or  periodical  type  (the  latter  type  is  valid  for  rectangular  domains  only).  Since  the  Laplace 
operator  with  the  mentioned  boundary  conditions  allows  only  a  discrete  spectrum,  the  traditional 
convolution  integral  of  the  Fourier  transform  should  be  replaced  with  summing  over  a  discrete  index. 
More  precisely,  for  a  complex-valued  function  A  =  A(x)  we  denote  its  direct  discrete  Fourier 
transfonn 


F {A}  I*?}  (flq, &2  vj ’•••)  ’ 


an  = 


-1 

Z2(Q) 


J  A(x) ■  en (x)dx , 


Q 


as  a  sequence  of  Fourier  coefficients  regarding  to  eigenfunctions  {en  (x)}^|  of  the  Laplace  operator 
meeting  the  corresponding  boundary  conditions:  Aj_e/7(x)  +  Xnen(x)  =  0 .  Similarly,  F  means  the 


+00 

inverse  discrete  Fourier  transform,  which  calculates  the  sum  F-1{a}=  Ya/^x)  .  In  such  a  way,  we 

n= 1 

denote  the  discrete  transfer  function  T  =  (pi, pn,...)  as  a  sequence  of  complex-valued 
multipliers,  which  directly  affect  the  corresponding  components  of  the  Fourier  spectrum.  In  our  terms, 
the  procedure  of  the  discrete  Fourier  filtering  of  the  function  A(x)  consists  of  calculating  its  discrete 

Fourier  transform  {a)=F{A},  component-wise  multiplying  with  the  transfer  function, 
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T  -a  =  (p\a\, P2CI2,—, pnan,...)  ,  followed  by  taking  the  inverse  discrete  Fourier  transform.  As  a  result, 
we  come  to  the  definition  of  discrete  Fourier  filtering  in  the  fonn  of  the  superposition 
Fj  {A}  =  F_1  {  T  •  F{A }  ]  that  is  placed  into  the  right-hand  side  of  equation  (2.1). 


Real-valued  parameters  8X  ,82,83  contained  in  the  right-hand  side  of  Eq.  (2.1)  are  responsible 

for  the  specific  arrangement  of  the  nonlinear  optical  system.  Namely,  equation  (2.1)  models  the 
dynamics  of  optical  systems  both  with  summation  of  the  input  and  the  feedback  light  field  and  with 
separation  of  these  fields.  In  the  former  case,  the  equation  allows  either  taking  into  account  the 
interference  of  these  fields  or  neglecting  it.  Moreover,  the  choice  of  these  parameters’  values  defines 
the  type  of  the  considered  nonlinear  medium  (focusing  or  defocusing).  For  example,  the  choice 
8}  =  S2  =  83  =  1  allows  us  to  simulate  an  optical  system  with  coherent  interaction  of  the  input  and 

feedback  fields  (the  interference  “is  on”),  while  the  choice  8X=8 3=1,  S2=  0  allows  us  to  neglect  the 

interference-type  interactions.  Should  we  set  <5,  =  S2  =0,  S3  =- 1 ,  then  (1)  will  simulate  an  optical  system 

with  separated  input  and  feedback  fields  in  the  case  of  a  defocusing  nonlinear  medium.  Feedback 
factor  K  is  proportional  to  the  intensity  of  the  input  light  wave. 

The  first  step  in  the  mathematical  study  of  the  optimization  problem  is  to  prove  the  existence  of 
the  solution  of  state  equation  (2.1)  in  the  appropriate  functional  space  for  all  the  admissible  filters.  To 
formulate  the  existence  and  uniqueness  theorem,  we  denote  by  a  space  of  complex-valued 

bounded  sequences  T  =  (pi,p2, ...,pn,...)  with  the  finite  norm  ||r||  =  sup  pn  | .  Using  f as  the 

n= in¬ 
admissible  set  allows  one  to  consider  a  wide  range  of  Fourier  filters  with  no  restrictions  of  smallness  to 
the  amplitudes  of  high  spatial  frequency  multipliers  to  be  supposed. 


1) 


Theorem  2.1.  Let  the  initial  condition  uq  e  LZ(D.)  and  the  input  light  wave 
Ain  e  H  ]  (Q)  n  C(Q) .  Then: 

For  any  discrete  filter,  T  e  £rrj,  on  any  finite  time  interval,  0  <t  <t\,  the  initial-boundary  value 
problem  for  equation  (2.1)  with  homogeneous  Dirichlet  or  Neumann-type  boundary  conditions  or 
periodic  boundary  conditions  has  the  unique  energy  class  solution 
u(x,t)  eL  (0,tp,H  (Q))nC([0,q];Z  (Q)).  The  solution  satisfies  the  initial  condition  in  the 
sense  of  traces,  the  equation  (2.1)  is  realized  in  the  sense  of  distributions. 

The  solution  obeys  the  estimation  \\u  ||^0  (  +  II w 


2) 


3) 


4) 


with 


|M|II2(0 , q  ;//*(□)) = 


the 

h 

(j(IK0 

0 


norms 


I V2  (0,q ; Hl  (Q))  -  C(t]  X1+  1 1  u°  I V2  (Q) 

11 U  1 1(0 -,L2 (Q)) =  ] 11  U{t)  1 V (Q)  ’ 


P2  + 1|  V±n(t)  ||22  )dt)^2 ,  where  the  constant  C(fi)>0  is 

L  (Q)  L  (Q) 


independent  of  u0 . 

The  solution  Lipschitz-continuously  depends  on  the  initial  function  and  the  Fourier  filter,  and  the 
estimation 

l|M""lf(0,<1;L2(n))+l!"“‘'llI2(0,(1;H1(Q))-C2<lll,0““0lk(Q)+l|r":r||f«.) 

holds  for  the  corresponding  solutions  u  =  w(x,  t\  Uq,T),  u  =  u(\,  t ;  Uq,T)  . 

Every  bounded  set  of  initial  data  is  attracted  to  some  set  M  (global  attractor),  which  is  compact 
in  L2  (Q)  and  has  a  finite  Hausdorff  dimension. 


2.2.2.  Statement  and  solvability  of  the  problem  of  optimal  Fourier  filtering. 
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The  optimization  problem  has  been  formulated  as  follows:  let  the  Fourier  filter  take  the  form 
Tqo  +  T  ,  where  T,tj  is  the  constant  part  of  the  filter  (for  instance,  =  0  ),  and  T  appears  as  a  variable 
one.  The  question  is  how  the  choice  of  the  specific  discrete  Fourier  filter  T  in  equation  (2.1)  allows 
one  to  obtain  at  a  given  point  of  time,  t  =  t\,  a  solution,  u(x,t  =  q ;  uq  ,  Trrj  +  T) ,  close  to  the  desired 

state,  rq(x).  In  the  notation  u(x,C,uq,  Trfj  +  T)  for  the  solution  of  problem  (2.1),  we  underline  its 
implicit  dependence  on  the  initial  phase  modulation  w0(x)  and  the  Fourier  filter  Tx  +  T  .  In  the  studies 
of  this  optimization  problem,  some  new  developments  have  been  achieved. 

Let  us  consider  a  linear  bounded  operator  D :  L1  (Cl)  — >  L1  (Q)  and  formulate  a  new 
minimization  problem  for  a  family  of  terminal  functionals 

J(T)  =  J  |  D(u(x,  t  =  ti;u0,Tao+T)~ui  (x))  |  dx  ->  inf,  feH,  (2.2) 

Q 

with  the  admissible  set  S  ,  which  takes  into  account  the  technical  restrictions  for  the  filters  considered. 
In  compare  with  the  “classical”  terminal  target  functional 

J\  (T)  =  J  |  u(x,t  =  t] ; uq ,  TrX}  +  T)  -  u\(x)  \ 2  dx  ,  the  operator  D  is  responsible  for  additional  features  of 
Q 

approximation  to  the  target  state  for  each  point  of  domain  Q  .  In  particular,  by  choosing  the  following 
operator 

D(v)  =  Jn,(x)(v(x)  -  |n'|_I  jnZn-(y)  ■  v(y)dy), 

we  can  embrace  a  practically  important  case  of  the  numerically  studied  problem  of  the  proximity  of  the 
sum  of  the  input  wave  phase,  <p(x),  and  the  additional  nonlinear  phase  modulation,  u(x,t) ,  to  the 

given  phase  distribution,  u\(x)  ,  in  the  bounds  of  the  controlled  subdomain  Cl'  at  the  time  point  t  =  : 

J  (T)  = 

=  j|zo'W(«(x;h;' Too  +  T)  +  (p(x)  -  M,(x)  -  |Qf jn-(y)(w(y;q;r00  +  T)  +  (p(y)  -  w^ylVyj  dx 
Q 

f  1,  if  xefi',  .  . 

where  Xn'(x)  =  \  ,  IQ  I  being  the  area  of  the  subdomain  Cl' . 

[0,  if  x  g  Q 

It  is  important  to  describe  possible  admissible  sets  S  for  which  the  optimization  problem  (2.1)- 

(2.2)  is  solvable,  i.e.  at  least  one  optimal  filter  T  satisfying  the  equality  J(T  )  =  inf  J(T )  does 

TeZ 

exist.  The  following  assertions  on  the  solvability  of  the  optimization  task  have  been  proven. 

Theorem  2.2.  Let  Ajn  e  C 1  ( Q ) ,  the  operator  D  be  a  linear  and  continuous  operator  in  Lr  (Cl), 
and  the  constant  part  of  the  Fourier  filter  Tra  e  be  fixed. 

1)  If  S  =  Hr,  where  HR  =  {T  =  (px,p2,..)  e  l2  : 1  Pk  1^  rkk\°^rk  <  R,k  =  1,2,..., R  >  0}, 

* 

then,  the  optimal  filter  set  HR  for  problem  (2.1)-(2.2)  is  not  empty, 

Hr  =  {T  e  HR  :  J(T  )  =  inf  J(T)}  ^  0  .  Moreover,  the  optimal  filter  set  HR  is  compact,  and  any 

TeHR 

sequence  minimizing  terminal  functional  (2.2)  converges  to  this  set. 

GO 

2)  If  S  =  Br,  where  BR-{T-  (p\, ,...)  e  f2  :  II  T  ||^  =  V|  pn  |2  <  R “}  ,  then,  the  optimal 

i2 

n= 1 

filter  set  BR  for  problem  (2. 1  )-(2.2)  is  not  empty,  BR  =  {T*  e  BR  :  J(T*)  =  inf  J(T)}  *  0  . 

TgBr 

* 

Moreover,  the  optimal  filter  set  BR  is  weakly  compact  in  the  space  f  2 . 

3)  For  the  gradient  J'(T)  ,  the  following  formula  is  valid: 
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I 

(J'(T),AT)  =  r~l  J  J (®r(«(x,0,7,oo  +  T),AT)£  y/(x,t)dxdt . 


where  (Ot,  AT)  =  2 K  B.Q\(S2Ain  +  S3ft  'Ain  exp(m)})  ■  FAT\Ain  exp (iu 

solution  of  the  conjugate  problem 

dy/  ?  * 

-r  —  +  tf/-ldA±i//  =  (Ou,^),  fe[0,fi], 
ot 

y/(x,  t  =  tx)  =  2D*  D(u(x,  t  =  ti;u0,Taa  +T)~  u\  (x)), 


and  y/{r,t )  is  the 


(; u , T),i//)  =  2KS2  ReU/„  •  i  •  exp(iu)F=  \Ainy/\  + 


+  2KS3  Re|A/w  •  i  •  exp (iu)Ff(FT  {Ain  exp (iu)}y/)j 

The  boundary  conditions  coincide  with  those  of  the  direct  problem  (2. 1). 

4)  Within  any  finite  time  interval,  0<  t<  t\,  initial-boundary  value  problem  (2.4)-(2.5)  with 

homogeneous  boundary  conditions  of  the  Dirichlet,  Neumann,  or  periodical  type  has  the  unique 
2  1  2 

solution  y/  e  L"(Q,t\,H  (Cl))  n  C([0,tj];Z  (Q)) ,  and  this  solution  obeys  the  estimation 

II  ^  Hc(0,/,;I2(O))  +  II  ^  llz.2(0,r1;ff1(O))~  F  H  ~  ^  Hi2(n)^’ 

where  the  constant  C  >  0  is  independent  of  ul . 


2.2.4.  Projection  finite-element  approximations 

To  solve  direct  problem  (2.1)  with  periodic  boundary  conditions  set  in  a  rectangle,  a 
rectangular  Cartesian  coordinate  system  is  usually  chosen.  An  implicit  projection-difference  scheme 
based  on  the  finite  element  method  has  been  proposed  and  developed.  We  used  space 

5*1  =  Lin{cpj  ,y(x)}  of  piecewise  continuous  linear  functions  (first-order  B-splines) 
<Pij(x)  =  (Pij(x \,x2)  =  <Pj(x j)  •  <Pj(x 2)  built  in  accordance  with  the  formulas 

<Pn  ixj)  =  Vo ( xj  ~  nhj),  <pN(xj )  =  <pQ (x-)  +  cpQ (xj  - 2k)  , 

(Pq(xj)  =  max{o,  1-  |  x ,  |  AJ1}  ,  n  =  1  -1,  j  =  1,2  . 

Then,  an  approximate  solution  was  being  sought  for,  as  the  element  w  =  (w°,...,wM)  e  S  of  the  space 
S  =  5)  x  5j  x  5)  x ...  x  S\  satisfying  the  projection-difference  scheme: 

M+ 1  times 

( T~l(wm  -  wm~l ),  <pj  +  lwm , <pj  + 1} (v j_wm ,  =  (®(wm ), <pj,  m  =  ,  V <p  e  Sj ,  (2.6) 

where,  as  before,  the  operator  <E>  stays  for  the  right-hand  side  of  (2.1).  The  angular  brackets  in  (2.6) 

9 

denote  inner  products  in  the  Hilbert  space  L  (Cl) .  In  the  software  developed,  we  used  a  more 
convenient  matrix  form  of  Eq.  (2.6): 

{(\  +  T)B  +  TDA)wm  =  T®(wm)  +  Bwm~l ,  m  =  (2.7) 

2 

where  A  =  B2AX  +  Bx  A2 ,  B  =  BXB2 ,  Bj  =  E  -  h  A/  /6 ,  E  is  the  unit  matrix,  and 

(A,/),  — i(/», -2 =  /, .  /», 0.....N,  - 1 ,  t  =  1,2 . 
n. 

The  above  scheme  (2.7)  is  implicit.  Therefore,  it  requires  iterations  in  order  to  obtain  the  solution  at  a 
new  temporal  layer.  Usually,  3-4  iterations  per  each  temporal  layer  are  enough  to  achieve  the  required 
accuracy. 

The  same  approach  is  applied  while  constructing  approximations  to  conjugate  problem  (2.4). 
Then,  an  approximate  solution  of  the  problem  (2.4)  is  being  sought  for,  as  the  element 
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w  =  )  e  S  satisfying  a  projection-difference  scheme  analogous  to  (2.6).  The  corresponding 

matrix  notation  looks  as  follows: 

((1  +  t)B  +  t  DA)wm -1  =  r  G{wm~l )  +  Bwm  ,  m  =  M,...,l .  (2.8) 

Scheme  (2.8)  is  being  solved  along  a  decreasing  sequence  of  indices.  Therefore,  this  scheme  is 
implicit.  The  choice  of  an  implicit  scheme  is  due  to  its  greater  stability  margin  with  respect  to  explicit 
schemes. 


2.2.5.  Tikhonov  regularization  method 


A  number  of  computer  experiments  were  performed  with  the  aim  of  minimizing  target 
functionals  of  the  considered  class  by  means  of  the  conditional  gradient  method.  On  this  way,  we 
noticed  that  after  a  certain  number  of  steps  further  iterations  did  not  considerably  decrease  the 
functional  value,  while  the  structure  of  the  optimal  filter  became  more  and  more  complex  (due  to  a 
larger  number  of  high-frequency  harmonics  to  transmit).  Such  an  unjustified  complication  of  the  filter 
structure  is  not  only  inconvenient  in  practice  but  also  indicates  certain  nonrobustness  of  the 
optimization  process.  To  make  the  optimization  robust,  we  used  Tikhonov's  regularization,  which  had 
been  earlier  applied  by  the  Team  members  to  abstract  extremum  problems  (A.N.  Tikhonov  and  F.P. 
Vasiliev  //  Banach  Center  Publications.  Vol.  3.  Mathematical  models  and  numerical  methods. 
Warszawa,  1978,  pp.  297-342;  F.P.  Vasiliev  and  M.  Yachimovich  //  Soviet  Mathematics  Doklady, 
Vol.  250,  No  2,  1980). 

We  have  implemented  a  version  of  the  Tikhonov  regularization  method  for  the  problem  of 
minimization  of  the  functional  J(T) ,  which  implies  constructing  a  family  of  Tikhonov  functionals 

Jk{T)  =  J{T)  +  ak\Tli, 


dependent  on  the  regularization  parameters  {a^  }^j  .  Then,  on  the  basis  of  the  expressions  for  the 
gradients,  J\{T)  =  J'(T)  +  2a/(T ,  for  every  k  =  1,2,...  one  step  of  the  conditional  gradient  method 


was  made  for  the  functional  J ^(T)  in  the  set  HR  .  As  a  result,  a  sequence  of  optimized  filters,  {Tk}  , 

was  obtained.  In  the  theory  of  regularization  of  the  conditional  gradient  method,  the  convergence  of 
the  iterations  to  the  set  of  optimal  controls  can  be  proven  under  rather  restricting  conditions,  such  as 
convexity  of  the  target  functional,  presence  of  a  saddle  node  for  the  Lagrange  function,  etc.  For  the 
Fourier  filtering  models  we  consider,  these  conditions  are  most  likely  not  met.  Flowever,  this  fact  does 
not  prevent  us  from  studying  possibilities  of  practical  application  of  iterative  regularization  methods 
and  robust  construction  of  a  sequence  of  optimized  filters.  The  properties  of  such  a  sequence  strongly 
depend  on  the  choice  of  the  regularization  parameters.  These  properties  were  numerically  studied  for 
Tikhonov  functionals  of  the  two  following  types. 

In  the  first  case,  all  the  regularization  parameters  were  equal:  {ak \kX  =  a  =  const  >  0  .  As  a 


result,  we  proceed  to  so-called  static  regularization  that  implies  adding  to  the  target  functional  a 
stabilizing  increment,  kY|  7’||“  ,  independent  of  the  iteration  number,  which  prevents  the  filter  norm 

from  increasing  during  the  optimization  process. 

In  the  second  case,  we  a  priori  constructed  a  power-function  sequence  of  the  regularization 
parameters,  {ak  =  C0  (k  + l)  r> ,  C0  >  (),/?  >  0 ,  which  was  often  used  in  regularization  methods.  As  a 

result,  we  obtained  a  particular  case  of  the  regularized  method  of  conditional  gradient. 

In  both  cases,  the  key  question  was  how  many  harmonics  of  the  discrete  Fourier  filter  should  be 
controlled  in  order  to  achieve  the  desirable  result  with  acceptable  accuracy.  Also,  we  compared  our 
numerical  data  with  the  results  obtained  in  the  case  of  no  regularization,  that  is,  when  {ak )'k=]  =  0 . 


Taking  one  group  of  such  experiments: 

hr={T  =  (php2,...):\pk\<l,k  =  \,2,-,N;pk=0,k>N  +  \}, 


To0=(  0,...,0 


N -times 


as  an  example,  let  us  describe  the  general  regularities  revealed. 
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For  the  phase  inhomogeneity  of  the  input  light  beam  having  the  form  of  ^»(x)  =  0.5cos(9x2) , 
we  took  the  target  function  rq(x)  shown  in  Fig  2.1a.  Fig  2.1b  demonstrates  the  total  phase  profile, 
u(\,ti;T )  +  (p(\) ,  before  the  optimization  is  started. 


Fig  2.1a.  Target  phase 
profde  U](x) . 


Fig  2.1b.  Total  phase 
profde  u(\,  t\ ;  T)  +  cp(x) 
before  the  optimization. 


In  the  case  of  controlling  a  limited  number  of  the  filter  components,  which  are  inside  a  desired- 
radius  circle  on  the  Fourier  plane,  close  results  were  obtained  both  for  iterative  regularization  with 
parameters  C0  =  1 ,52588e-5,/f  =  0.45 (Fig  2.2)  and  in  the  experiment  with  no  regularization,  that  is, 

f°r  {a,  }ITi  =  0  (Fig  2.3). 

In  the  experiment  with  static  regularization  with  parameters  {ak =a  =  1.52588e-5,  the 

approximation  (see  Fig  2.4a)  appeared  to  be  less  successful  with  respect  to  those  shown  in  Fig  2.3a 
and  Fig  2.2a.  This  can  be  explained  by  the  fact  that  the  resultant  optimized  control  shown  in  Fig  2.4b 
has  a  smaller  nonn  as  compared  with  the  controls  obtained  in  the  previous  experiments  (Fig  2.2b  and 
Fig  2.3b). 

Similar  regularities  were  revealed  during  a  series  of  computer  experiments  where  all  the 
available  components  of  the  discrete  Fourier  filter  were  controlled.  In  these  experiments,  good 
approximations  to  the  target  function  ii\ (x)  were  obtained  both  with  iterative  regularization  with 

parameters  C0  =  1 .52588c - 6,/f  =  0.25  (Fig  2.5)  and  with  no  regularization,  that  is,  for  \ak  \'k=]  =  0  (Fig 

2.6).  Nevertheless,  these  experiments  qualitatively  differ  from  the  experiments  with  a  limited  number 
of  controlled  harmonics.  This  difference  is  in  the  fact  that  the  iterative  regularization  performed  for  all 
the  controlled  filter  components  allowed  us  to  clearly  demonstrate  the  possibility  of  obtaining  a  good 
approximation  with  simultaneously  decreasing  the  number  of  controlled  harmonics  (as  compared  with 
the  experiment,  in  which  no  regularization  was  performed).  In  this  case,  it  becomes  possible  not  to 
control  a  certain  number  of  harmonics  belonging  to  peripheral  portions  of  the  filter  (in  Fig  2.7,  these 
harmonics,  which  are  ignored  in  the  filter  shown  in  Fig  2.5b  with  respect  to  the  filter  shown  in  Fig 
2.6b,  are  red-highlighted).  This  fact  emphasizes  a  positive  effect  of  Tikhonov’s  iterative  regularization 
for  obtaining  both  efficient  and  non-complex  Fourier  filters. 
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Fig  2.2b.  Modulus  of  the 
discrete  filter  T  . 


Fig  2.2a.  Total  phase 
profde  u(x,  t\ ;  T)  +  qy{x) 
after  optimization  with 
iterative  regularization 


Fig  2.2c.  Argument  of  the 
discrete  filter  T . 


Fig  2.3a.  Total  phase 
profde  u(x,ti;T)  +  (p(x ) 
after  optimization 
without  regularization 


Fig  2.3b.  Modulus  of  the 
discrete  filter  T  . 


Fig  2.3c.  Argument  of  the 
discrete  filter  T  . 


Fig  2.4a.  Total  phase 
profde  u(x,  t\ ;  T)  +  <p(x ) 
after  optimization  with 
static  regularization 


Fig 


2.4b.  Modulus  of  the  Fig  2.4c.  Argument  of  the 

discrete  filter  T  .  discrete  filter  T  . 
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Fig  2.6a.  Total  phase 
profile 

u(x,t\,T)  +  cp(x )  after 
optimization  without 
regularization 


Fig  2.6c.  Argument  of  the 
discrete  filter  T  . 


Fig  2.7.  Red  points  mark  those 
harmonics  processed  by  the 
filter  shown  in  Fig  2.6b  and 
ignored  in  the  filter  shown  in 

Fig  2.5b 


2.2.6.  Dependence  of  the  optimization  quality  on  the  number  of  controlled  harmonics 


Some  new  statements  were  considered  for  the  problem  of  forming  a  given  phase  distribution  in 
a  light  wave  by  means  of  optical  feedback  systems  with  Fourier  filters,  which  correspond  to  model 
(2.1)  for  the  following  values  of  right-hand  side  parameters:  Sl  =  82  =0,  S3  =-l.  The  optimization 
quality  was  estimated  by  the  following  functional: 

Jk(T)  = 


u(^;tl;Tx  +  T)  +  (p(x)-ul(x) 


+  ak 


Jn  0(y)(«(y;  tl;Ta)+T)  +  <p(  y)  -  ux  (y))dy  ' 
la0(y)dy 


2 

dx  + 


In  the  numerical  simulation  series  we  performed,  the  following  parameters  of  iterative  regularization 
were  chosen:  C0  =  0.833/V  ',/f  =  0.45,  where  N  is  the  number  of  controlled  harmonics,  and  the 
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sequence  {ak  \k  \  was  constructed  according  to  the  formula  \ak  jy*  =  C0  (k  + 1)  13 .  The  goal  of  the 

numerical  simulations  was  to  figure  out  the  number  of  controlled  harmonics,  enough  to  achieve  a  good 
quality  of  forming  the  desired  phase  distribution.  The  key  parameters  were  the  following:  R ,  the 
radius  of  a  circle  centered  with  respect  to  the  Fourier  plane,  which  defined  the  number  of  "active" 
harmonics,  N{R )  (all  the  harmonics  being  inside  the  circle).  The  two  additional  quantities,  ?]((p)  and 

77(1^)  ,  we  used  to  explain  the  results  obtained,  were  the  fractions  of  the  power  spectra  of  (respectively) 
the  input  field  Ain (x)  =  dm(x)|exp{/7/9(x)j  and  the  target  field  \Ain  (x)|  £xp f/w,  (x) } ,  which  lay  within  the 
controlled  zone  (the  circle  of  the  radius  R).  The  control  quality  was  estimated  with  the  use  of  a  smooth 
weight  function,  #(x)  e  C°°  (Q) .  Taking  one  group  of  such  simulations  as  an  example,  let  us  describe 
the  regularities  revealed. 

For  the  phase  inhomogeneity  of  the  input  light  beam  having  the  form  of  (p(x)  =  0.5  cos(16x2)  , 
we  took  the  target  function  ux  (x)  =  exp(-  5((x,  -  n  f  +  (x2  -  n  f ))  shown  in  Fig  2.8a. 


Fig  2.8a.  Target  phase  Fig  2.8b.  Total  phase  profile 
profile  u\  (x) .  u(x,  q ;  T)  +  (p{x )  before  the 

optimization. 


For  those  simulations  where  the  controlled  zone  radius,  R ,  was  large  enough  to  fulfill  the 
conditions  r/((p)  >  0.99  and  //(u, )  >  0.99,  the  optimization  results  were  the  best  out  of  the  series.  In 
this  case,  the  number  of  controlled  harmonics  were,  of  course,  quite  large,  while  the  rest  portion  of  the 
Fourier  filter  was  just  opaque  (Fig  2.9).  ft  can  be  seen  from  Fig  2.9b  that  only  those  harmonics 
belonging  either  to  the  spectrum  of  the  input  wave  or  to  that  of  the  target  profile  take  part  in  the 
optimization  process.  This  feature  is  important,  because  it  may  help  with  adaptively  choosing  the 
structure  of  the  controlled  set  of  Fourier  filters  for  any  specific  input  data. 


Fig  2.9a.  Total  phase  profile 
u(x,  q ;  T)  +  (p{x )  after  the 
optimization.  R  —  50 . 


Fig  2.9b.  Modulus  of  the 
discrete  filter  T . 
N(R)  =  7845  . 


Fig  2.9c.  Argument  of  the 
discrete  filter  T  . 
N(R)  =  7845  . 
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When  decreasing  R  ,  we  revealed  the  following  regularity.  In  those  simulations  where  the  controlled 
zone  radius  were  still  large  enough  to  fulfdl  the  conditions  r/(<p)  >  0.92  -s-  0.93  and 
77(1^)  >  0.92  -5-  0.93 ,  the  optimization  result  remained  rather  good  (Fig  2.10),  though  distinctly 
conceding  to  those  simulations  where  rj{cp)  >  0.99  and  >  0.99  (Fig  2.9).  When  we  decreased  R 
in  such  a  way  that  to  fulfill  /]((p)  <  0.90  or  //(u, )  <  0.90 ,  the  optimization  quality  considerably 
degraded  with  respect  to  the  previous  simulation  series  (compare  Fig  2.11  and  Fig  2.9-Fig  2.10). 


Fig  2.10a.  Total  phase 
profde  u(x,  t\ ;  T)  +  cp{x) 
after  the  optimization. 
R  =  19. 


Fig  2.10b.  Modulus  of  the 
discrete  filter  T . 
N(R)  =  1129  . 


Fig  2.10c.  Argument  of  the 
discrete  filter  T  . 

N(R)  =  1 129  . 


Fig  2.11a.  Total  phase 
profde  u(x,ti,T)  +  (p(x) 
after  the  optimization. 

R  =  14. 


Fig  2.11b.  Modulus  of  the 
discrete  filter  T  . 
N(R)  =  613 . 


Fig  2.11c.  Argument  of  the 
discrete  filter  T  . 
N(R)  -  613 . 


Further  decreasing  R,  in  such  a  way  that  rj((p)  <  0.88  or  //(u, )  <  0.88 ,  we  observed  further 
deterioration  of  the  optimization  quality  (Fig  2.12)  with  respect  to  the  above  series  (Fig  2.9-Fig  2.11). 
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Fig  2.12a.  Total  phase 
profile  u(\,ti,T)  +  <p(x)  after 
the  optimization. 

R  =  6. 


Fig  2.12b.  Modulus  of  the 
discrete  filter  T  . 
N(R)  =  113. 


Fig  2.12c.  Argument  of  the 
discrete  filter  T  . 
N(R)  =  113. 


Fig  2.13  presents  the  dependence  of  the  optimized  value  of  the  target  functional  J  on  the 
radius  R  .  The  numerical  data  on  which  this  plot  is  based  are  listed  in  the  table  below  (Fig  2.14). 


Fig  2.13.  Target  functional  value  after  the  optimization 
versus  the  radius  of  the  controlled  zone. 
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R 

N(R) 

J 

77(Wi) 

1 

5 

0,7134 

47,71 

47,58 

2 

13 

0,6414 

75,76 

75,45 

3 

29 

0,4839 

82,94 

85,43 

4 

49 

0,3789 

84,52 

92,29 

6 

113 

0,3636 

87,25 

98,58 

9 

253 

0,3411 

87,99 

99,88 

12 

441 

0,3231 

88,12 

99,98 

14 

613 

0,2780 

88,50 

100,00 

15 

709 

0,2264 

89,95 

100,00 

16 

797 

0,1375 

93,43 

100,00 

17 

901 

0,0767 

97,37 

100,00 

19 

1129 

0,0530 

99,60 

100,00 

35 

3853 

0,0513 

100,00 

100,00 

50 

7845 

0,0482 

100,00 

100,00 

Fig  2.14 


The  main  conclusion  we  may  derive  from  the  above  results  is  the  following.  To  achieve  the  best 
optimization  quality,  it  is  necessary  to  choose  such  a  number  of  harmonics  to  control  that  to  fulfill  the 
conditions  rj{(p)  >  0.99  and  //(u, )  >  0.99 ,  meaning  that  the  controlled  zone  of  the  Fourier  plane  should 
contain  al  least  99%  of  the  power  of  both  the  input  and  the  target  field  spectrum. 

2.2.7.  Dependence  of  the  optimization  quality  on  the  parameter  K 

To  figure  out  how  the  parameter  K  (proportional  to  the  intensity  of  the  input  wave)  affects  the 
quality  of  forming  the  desired  phase  distribution  in  the  output  wave,  we  performed  some  numerical 
simulations,  taking  the  following  values  of  the  iterative  regularization  parameters: 
C0  =  0.833AT1,/?  =  0.45  ,  where  N  =  1793  is  the  number  of  controlled  harmonics,  and  the  sequence 

\ak  \k=\  was  constructed  according  to  the  formula  \ak  =  C0  (k  + 1)  r> .  The  optimization  problem 

was  solved  with  the  conditional  gradient  method,  in  which  a  combined  stopping  criterion  was  used  to 
take  into  account  both  the  discrepancy  value  and  the  number  of  executed  iterations  of  the  gradient 
method.  The  control  quality  was  estimated  using  a  smooth  weight  function,  6(\)  e  C“(Q).  Thus,  all 
the  numerical  simulation  results  shown  below  were  obtained  at  the  same  fixed  parameters  of  the 
optimization  method.  Let  us  describe  the  regularities  revealed,  using  one  series  of  simulations  as  an 
example. 

In  each  numerical  experiment  of  this  series,  the  diffusion  coefficient  was  set  to  l2d  =  0.0025  ,  while 
the  parameters  of  the  input  field  A  m  ( x  )  =  \Ajn  ( x  ) I  exp {  i<p  (  x)\  were  the  following: 
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[  1,  if  xe  Q', 

A.  (x)  e  C°°  (Q) ,  4„(x)  =  i  where  Q'cfi1',  (p(x)  =  2cos(7x2  -  3Xj)  .  The  target  profile 

[0,  if  x  £  Q  , 

was  defined  by  the  formula  u t(x)  =  2cos(2;q  +  3x,  )cos(4x2  -5jCj)  (see  Fig  2.15a.): 


Fig  2.15a.  Target  phase 
profde  U\  (x) . 

Fig  2.15b.  Total  phase 
profde  u(x,ti,T)  +  (p{x) 
before  the  optimization. 

A  =  0.5. 

Fig  2.15b  shows  the  total  phase  profile  u{x,tx\T )  +  <p(x)  before  the  optimization,  for  K  -  0.5  . 
We  should  stress  that  the  2-radian  amplitude  of  both  the  input  phase  distortion  and  the  target  phase 
profile  is  quite  large.  Fig  2.16  clearly  demonstrates  that  K  =  0.5  is  absolutely  not  enough  to  achieve 
good  optimization  quality.  Indeed,  at  this  value  of  K  no  considerable  approach  to  the  target  profile  is 
obtained  (compare  Fig  2.15a.  and  Fig  2.16a.),  and  the  final  result  of  the  optimization  procedure 
weakly  differs  from  the  initial  phase  profile  (Fig  2.15b).  Fig  2.16b  shows  that  almost  all  the  harmonics 
are  equally  (non-adaptively)  transmitted  by  the  filter.  This  fact  clearly  indicates  that  at  K  =  0.5  the 
depth  of  the  additional  phase  modulation  u(x,t{\T)  is  not  enough  to  efficiently  suppress  the  distortion 
and  simultaneously  form  the  target  profile  shown  in  Fig  2.15a. 


Fig  2.16a.  Total  phase 
profde  u(x,  t\ ;  T)  +  (p(x) 
after  the  optimization. 
K  =  0.5 


Fig  2.16b.  Modulus  of  the 
discrete  filter  T . 

K  =  0.5 . 


Fig  2.16c.  Argument  of  the 
discrete  filter  T  . 

A  =  0.5. 


In  the  next  numerical  experiment,  the  input  wave  intensity  was  increased  to  K  =  1.1.  The 
results  of  this  experiment,  shown  in  Fig  2.17,  demonstrate  that  the  optimization  quality  considerably 
improved.  It  can  be  seen  that  the  fdter  obtained  for  K  =  1.1  (see  Fig  2.17b,c)  is  less  trivial  and  more 
adaptive  than  the  previous  one  (Fig  2.16b,c).  The  latter  feature  is  confirmed  by  the  fact  that  the  total 
phase  profile  after  the  optimization  (Fig  2.17a)  decently  fits  the  target  profile  (Fig  2.15a). 
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Fig  2.17b.  Modulus  of  the 
discrete  filter  T . 

K  =  1.1. 


Fig  2.17c.  Argument  of  the 
discrete  filter  T  . 

K  =  \.\. 


Fig  2.17a.  Total  phase 
profde  u(x,  t\ ;  T)  +  qj{x) 
after  the  optimization. 
K  =  1.1. 


Similar  regularities  were  revealed  for  K  -  3.2  (Fig  2.18)  and  A  =  7.0  (Fig  2.19). 


Fig  2.18a.  Total  phase 
profde  u(x,  t\ ;  T)  +  qj{x) 
after  the  optimization. 


Fig  2.18b.  Modulus  of  the 
discrete  filter  T . 

K  -3.2. 


Fig  2.18c.  Argument  of  the 
discrete  filter  T  . 

A  =  3.2. 


K -3.2. 


Fig  2.19a.  Total  phase 
profde  u(x,  t\ ;  T)  +  cp{x) 
after  the  optimization. 
A  =  7.0. 


Fig  2.19b.  Modulus  of  the 
discrete  filter  T . 


Fig  2.19c.  Argument  of  the 
discrete  filter  T  . 


A  =  7.0. 


A  =  7.0. 


Fig  2.20  shows  a  plot  of  the  target  functional  values  obtained  as  a  result  of  the  optimization 
versus  the  parameter  K .  The  plot  illustrates  the  following  feature.  Choosing  larger  values  of  K  allows 
one  to  synthesize  a  Fourier  filter  that  is  better  adapted  to  the  specific  parameters  of  the  input  and  the 
target  field.  For  larger  values  of  K,  the  optimization  results  are,  as  a  rule,  not  worse  than  those  for 
smaller  values  of  the  parameter.  This  feature  is  due  to  the  fact  that  a  larger  value  of  K  leads  to  a  larger 
depth  of  the  additional  phase  modulation  ufxJ^T)  and,  in  turn,  to  wider  possibilities  of  correcting  the 
input  phase  profde  (p{x)  with  simultaneously  forming  the  desired  phase  profde  iq(x) . 


K 


Fig  2.20.  Target  functional  value  versus  the  parameter  K. 

Thus,  to  achieve  better  optimization  quality,  it  is  advisable  to  choose  a  rather  large  value  of  the 
feedback  factor  K . 

2.2.8.  Dependence  of  the  optimization  quality  on  the  diffusion  coefficient  of  the  nonlinear 

medium  particles. 

To  figure  out  how  the  effective  diffusion  coefficient,  D  =l2d,  of  the  nonlinear  medium  particles 
affects  the  quality  of  forming  the  desired  phase  distribution  in  the  output  wave,  we  performed  some 
numerical  simulations,  taking  the  following  values  of  the  iterative  regularization  parameters: 
C0  =  0.8337V1,/?  =  0.45  ,  where  N  =  4293  is  the  number  of  controlled  harmonics,  and  the  sequence 

\ak  }+”  was  constructed  according  to  the  formula  {ak  =  C0 {k  +  \)  /; .  The  optimization  problem 

was  solved  with  the  conditional  gradient  method,  in  which  a  combined  stopping  criterion  was  used  to 
take  into  account  both  the  discrepancy  value  and  the  number  of  executed  iterations  of  the  gradient 
method.  The  control  quality  was  estimated  using  a  smooth  weight  function,  0(x)  e  C°°(Q) .  Thus,  all 
the  numerical  simulation  results  shown  below  were  obtained  at  the  same  fixed  parameters  of  the 
optimization  method.  Let  us  describe  the  regularities  revealed,  using  one  series  of  simulations  as  an 
example. 

In  each  numerical  experiment  of  this  series,  the  feedback  factor  was  set  to  K  =  3,2 ,  while  the 
parameters  of  the  input  field  Ain  (x)  =  |  A  jn  ( x )  |  exp {  i<p  (x  )\  were  the  following: 
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K (x)  e  c°° (Q) >  An (x)  =  ' ’  ^  X  ^  where  Cl'  <=  Cl" ,  <p(x)  =  2cos(8(Xl  +  x2 ))cos(8(x,  -x2)).  The 

[0,  if  x?n , 

target  profile  is  shown  in  Fig  2.21a. 


Fig  2.21a.  Target  phase  Fig  2.21b.  Total  phase 

profile u\ (x) .  profile  u(x,ti;T)  +  <p(x ) 

before  the  optimization. 
D  =  0,0055 . 


Fig  2.21b.  shows  the  total  phase  profile  u(x,t{;T )  +  <p(x)  before  the  optimization,  for  K  -  0.5  . 
We  should  stress  that  the  2-radian  amplitude  of  both  the  input  phase  distortion  and  the  target  phase 
profile  is  quite  large. 

For  the  first  numerical  experiment  of  this  series,  we  chose  D  =  0,0025 .  The  results  shown  in 
Fig  2.22  demonstrate  rather  good  quality  of  forming  the  desired  phase  profile  as  well  as  a  complex 
adaptive  structure  of  the  optimal  filter  (Fig  2.22b, c). 


Fig  2.22a.  Total  phase 
profile  u(x,  t] ;  T)  +  cp(x) 
after  the  optimization. 
D  =  0,0025 . 


Fig  2.22b.  Modulus  of  the 
discrete  filter  T 
D  =  0,0025 . 


Fig  2.22c.  Argument  of  the 
discrete  filter  T  . 

D  =  0,0025 . 


Increasing  the  diffusion  coefficient  D  considerably  decreases  the  quality  of  the  numerical 
optimization.  This  is  caused  by  strengthening  of  diffusive  coupling  between  the  nonlinear  medium 
particles.  At  the  same  time,  the  larger  the  diffusive  coefficient,  the  higher  the  average  transmittance  of 
the  optimal  filter  (see  Fig  2.22b,  2.23b,  2.24b,  2.25b). 
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Fig  2.23a.  Total  phase 
profile  u(x,  t\ ;  T)  +  cp{x) 
after  the  optimization. 
D  =  0,0040 . 


Fig  2.23b.  Modulus  of  the 
discrete  filter  T . 

D  =  0,0040 . 


Fig  2.23c.  Argument  of  the 
discrete  filter  T  . 

D  =  0,0040 . 
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Fig  2.26  shows  a  plot  of  the  target  functional  values  obtained  as  a  result  of  the  optimization 
versus  the  diffusion  coefficient  D .  The  plot  demonstrates  the  following  regularity:  increasing  the 
diffusion  coefficient  sufficiently  impairs  the  optimization  quality.  Since  the  diffusion  coefficient 
defines  the  transverse  resolution  of  the  nonlinear  medium,  it  can  be  concluded  that  in  practice  it  is 
advisable  to  use  media  with  a  high  spatial  resolution.  For  example,  it  can  be  a  high-resolution  optically 
addressed  liquid  crystal  light  valve. 


D 

Fig  2.26.  Target  functional  value  versus  the  diffusion  coefficient. 


2.2.9.  Study  of  the  optimization  quality  achieved  in  optical  systems  with  two  feedback  loops. 


In  this  Section,  we  consider  a  theoretical  model  of  an  optical  feedback  system,  in  which  the 
additional  phase  modulation  u  =  u{x,y,t )  of  the  light  wave  passed  through  the  nonlinear  medium  is 
governed  by  the  equation 

du  1 

t - \-u  —  DA±u  =  I  fo  +  c\l  fbdt' ,  I  fh  =  K[\  +  y  cos(u  +  (p)],  (2.9) 

dt  o 

Here,  T  is  the  response  time  of  the  nonlinear  medium,  D  is  the  diffusion  coefficient  of  its  particles, 
K  is  the  feedback  factor,  dependent  on  the  nonlinear  properties  of  the  medium  and  the  intensity  of  the 
input  wave,  y  is  a  coefficient  dependent  on  the  average  intensity  of  the  feedback  wave,  (p  is  the  phase 

profile  of  the  input  wave,  and  C  is  the  factor  of  integral  feedback.  Equation  (2.9)  is  to  be  supplemented 
with  appropriate  initial  and  boundary  conditions. 

At  c  =  0 ,  the  above  model  describes  the  dynamics  of  a  number  of  well-known  nonlinear 
optical  systems,  such  as  a  nonlinear  interferometer  with  optical  feedback  or  a  nonlinear  single-pass 
ring  resonator.  It  is  known  that  such  systems  are  capable  of  automatically  suppressing  the  input  phase 
distortions  (p .  In  addition,  active  control  of  the  light  wave  phase  by  means  of  a  digital  processor 

connected  to  the  feedback  loop  can  considerably  improve  the  suppression  quality. 

At  the  same  time,  it  is  more  practical  and  effective  to  use  purely  optical  control  provided  by 
nonlinear  media.  It  is  supposed  that  it  is  possible  to  improve  the  quality  of  phase  distortion  suppression 
(with  respect  to  the  case  c  —  0 )  by  supplementing  the  optical  system  with  an  additional  feedback  loop, 
whose  functionality  corresponds  to  the  integral  term  in  the  right-hand  size  of  Eq.(2.9).  Our  numerical 
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results  confirm  this  supposition.  For  instance,  Fig.2.27  shows  the  total  phase  profile,  u  +  (p ,  at  t=  0  for 
the  case  when  (p  is  a  2D  phase  grating  with  an  amplitude  of  about  2  rad.  Fig.2.28  demonstrates  the 

effect  of  suppression  of  these  phase  distortions  for  c  =  0 .  It  can  be  seen  from  Fig.2.28  that,  as  time 
goes  on  ( t  =  2r ),  the  phase  profile  u  +  (p  gets  somewhat  flattened.  However,  Fig.2.28  clearly 

indicates  that,  under  the  same  conditions,  the  integral  feedback  (c  =  1)  allows  us  to  achieve  much 
better  results.  Moreover,  this  feature  is  observed  for  a  wide  range  of  system's  parameters.  Solving 
Eq.(2.9)  over  a  longer  time  interval,  one  could  see  that  the  total  phase  profile  u  +  (p  gradually  tends  to 

a  stationary  one  and  simultaneously  continues  to  flatten.  Such  a  dynamics  demonstrates  that  the  system 
can  be  applied  to  phase  distortion  suppression  tasks.  From  this  point  of  view,  it  is  necessary  to  further 
investigate  in  detail  how  the  suppression  quality  depends  on  system's  parameters.  It  would  be  also 
useful  to  consider  the  case  of  non- stationary  phase  distortions  . 


Fig.  2.27.  Total  phase  profile  U  +  (p  at  the 
moment  t=  0. 


Fig.  2.28.  Total  phase  profile  U  +  (p  at  the 
moment  t  —  'll  (c  =  1 ). 


Fig.  2.28.  Total  phase  profile  U  +  (p  at  the 
moment  t  =  2r  (C  =  0 ). 
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2.3.  Controlling  Irreversible  Transforms  of  Spatial  Arguments 

The  statement  of  optimal  argument  transformation  g(x)  =  (gi(x),g2(x)) :  £2  — »  R  has  been 
proposed  for  the  model  of  2D  feedback  optical  systems,  governed  by  the  functional-differential 
diffusion  equation 

3u 

—  +  u  -  DA±u  =  Fg(u),  te[0,T],  u(x,t  =  0)  =  w0(x),  «  |5nx(0,n=  0  •  (3-1) 

Here,  x  =  (x\,X2)  is  the  transverse  spatial  argument  varying  in  spatial  domain  Q  with  smooth 
boundary  dO.  or  in  a  rectangle;  u  =  u(x,t)  =  u(x,t;  g)  is  the  additional  phase  modulation  in  the 
nonlinear  medium  layer;  is  the  transverse  Laplace  operator;  D  is  the  diffusion  coefficient;  f  g  ( u ) 

is  the  superposition  of  the  nonlinear  term  F(u )  and  a  controllable  transform,  g,  of  the  spatial 
argument.  For  instance,  the  formula  F(w)  =  K(\  +  y  cos  u )  is  used  for  modeling  field  interactions  in 
an  interferometer  or  passive  ring  resonator  with  a  thin  Kerr-type  medium,  factor  K  is  proportional  to 
the  intensity  of  the  input  light  wave,  y  is  the  visibility  of  the  interference  pattern. 


2.3.1.  Distributed  statement  of  irreversible  transform  of  spatial  argument 


In  the  case  of  a  reversible  transform,  g(x)  =  (g1(x),g2(x))  :  Q — onto  >Q  ci  R1 ,  the  classical 
formula 

Fg(u)(x,t)  =  F  O(g-1(x),0)  (3.2) 

is  used.  It  describes  a  point-to-point  superposition  of  F (u)  with  the  inverse  transform  g  1  ( x )  .  For 
instance,  reflection  and/or  rotations  in  a  circular  domain  are  valid  transforms  that  are  easy  to  realize  by 
means  of  reflecting  mirrors  and  rotating  prisms. 

A  complicated  mathematical  problem  is  how  to  give  an  adequate  interpretation  of  such  a 
superposition  for  an  arbitrary  Lebesgue-measurable  transform.  We  propose  to  consider  such  a 
superposition  in  the  generalized  sense,  i.e.  as  a  distribution  defined  by  the  formula 


(^g(M),^}  =  |  F{u{x,t))  (p(g(x))  dx  ,  V^eC(Q).  (3.3) 

The  advantage  of  this  approach  is  that  there  is  no  inverse  transform  implemented  in  the  integral 
calculation  and  only  a  direct  transform  acts  in  the  superposition  with  the  test  function  cp{x).  It  is 
important  that  for  the  class  of  reversible  transforms,  the  distribution  in  (3.3)  is  regular  and  coincides 
with  that  of  (3.2).  The  following  new  result  on  solvability  and  smoothness  has  been  proved. 

Theorem  3.1.  Let  the  initial  condition  w0eZ2( Q)  and  the  function  F  be  bounded  and 
continuously  differentiable  with  a  bounded  derivative.  Then,  for  any  Lebesgue-measurable  transform 
geG,  the  conjugate  problem  (3.1)-(3.3)  has  the  unique  generalized  solution 

u  e  L2  (0,  T;  Hi  s  (Q)),  dtu  e  L~  (0,  T;  H  S(Q))  with  arbitrary  s  e  (1,3/2).  The  solution  Holder- 


continuously  depends  on  the  argument  transforms  in  the  sense  of  the  estimate 


\u  -  u 


h\o,T-H2-s(m+^u~u 


1 Lp{QT ) 


^c(p)\\g-g 


II P 

NL2(Q) 


with  p  =  4p  1  -  1  and  arbitrary  p  e  (2,4) . 


2.3.2.  Statement  of  optimal  control  of  argument  transforms 

The  goal  of  the  optimization  is  to  approximate  to  a  given  function,  u\(x,t) ,  the  solution  of 
(3.1).  The  target  functional  estimates  the  quality  of  the  approximation  being  the  optimization  result. 
This  functional  is  determined  by  the  requirements  of  the  problem  formulated  and  usually  coincides 
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with  the  norm  of  deviation  from  the  assigned  target,  measured  in  a  certain  functional  space.  It  is 
proposed  to  consider  a  weighed  integral  functional  of  the  following  kind: 

J(g)  =  \  p(x,t)\u(x,t;g)-ui(x,t)\2  dxdt ,  (3.4) 

Qt 

which  allows  embracing  different  statements  of  optimization  problems.  In  such  a  way,  varying  the 
form  of  the  function  p(x,  t )  makes  it  possible  to  consider  the  solution  approximation  within  the  whole 
time  interval  and  spatial  domain,  as  well  as  on  its  certain  portion  and  subdomain.  Parameter  p(x,  t ) 
determines  the  relative  contribution  of  each  temporal  layer  in  the  above  functional.  The  minimization 
problem  for  the  functional  /( g)  is  fonnulated  on  the  admissible  set  G  of  Lebesgue-measurable 
transforms  satisfying  pointwise  restrictions: 

/(g)  — »  inf ,  G  =  {ge  Z2( Q) :  g(Q)  e  Q}  .  (3.5) 

G 

Theorem  3.2.  Let  the  conditions  of  theorem  3.1  hold,  the  target  function  «,  e/2( Q) ,  and 
0<  p  e  Lf  (Q) .  Then  the  target  functional  /(g)  is  correctly  defined  on  the  admissible  set  G  and 
Holder-continuously  depends  on  the  argument  transforms. 

2.3.3.  Conjugate  problem  and  target  functional  gradient 

To  compute  the  transform  sequence  that  minimizes  the  functional,  we  have  used  a  variation 
approach  based  on  extracting  the  main  linear  part  of  the  increment  of  functional  /(g) .  Let  <//  =  <//(x,/) 
denote  the  solution  of  so-called  conjugated  problem: 

-dti//  +  i//-DA1i//  =  2p(u-ul)  +  F\u)i//(g(x),t),  ^  lanx(o,n=  H=r=0-  (3-6) 


Note  that,  unlike  that  in  Eq.  (3.1),  in  the  right-hand  side  of  Eq.  (3.6)  the  sought  function  is  superposed 
to  transfonn  g(x)  in  the  classical  sense.  Furthermore,  transfonn  g(x)  is  implicitly  contained  in  (3.6) 
through  the  solution  u  =  u(x,t;  g)  of  the  diffusion  equation  (3.1). 

Theorem  3.3.  Let  the  conditions  of  theorems  3. 1-3.2  hold  and  the  target  function  have 

additional  property,  u\  e  Z4  ( Qj ) .  Then,  for  any  Lebesgue-measurable  transform  ge  G  the  conjugate 

2  1 

problem  (3.6)  has  the  unique  solution  i//  e  H p  (Qp)  with  arbitrary  p  e  (2,4) ,  and  this  solution  obeys 


the  estimation  C(p)  ,  where  the  constant  is  independent  of  g  .  The  solution  Holder- 

continuously  depends  on  the  argument  transforms  in  the  sense  of  the  estimate 


\\V'-ft\\H2,\QT)+\\V/-ft 


G([0,T]-Hlcim 


^C{p)\\g-g\f2 


'Zr(Q) 


with  (3  =  4p  1  -  1  and  arbitrary  q  e  (2,+oo) . 

Theorem  3.4.  Let  the  conditions  of  theorem  3.3  hold.  Then,  for  any  Lebesgue-measurable 
transform  geG,  the  target  functional  /(g)  is  uniformly  quasi-differentiable,  i.e.  the  following 
increment  formula  is  valid: 

/(g  +  h)-/(g)  =(j'(g),h)i2(Q)+i?(g,h),  Vg,g  +  heG, 
with  the  reminder  term  R{ g,  h)  ,  which  is  unifonnly  estimated  on  G  with  the  inequality 

|tf(g,h)|  <  C||h|f^n),//>0. 

The  gradient  J'(g)  e  L'  (Q,)  is  calculated  according  to  the  formula 
T 

J'(g)  =  {  F{u{x,  t;  g)} Vy  (g(x),  t)dt  (3.7) 

0 


37 


iil/3 


and  meets  the  Holder-continuity  property,  ||  J’(g)  -  J  ( g )  JJ^^.  ^  L  ||  g  —  g  n/2(.Q) . 


2.3.4.  Description  of  gradient-type  methods 

To  numerically  implement  the  method  based  on  the  previously  derived  gradient  fonnula  and  on 
the  additional  information  on  its  Holder-continuity,  three  methods  were  elaborated.  The  first  one,  a 
special  gradient  projection  method,  was  elaborated  and  theoretically  analyzed  in  the  framework  of 
the  Project.  In  this  version,  the  sequence  of  optimizing  transforms  is  derived  from  the  following 
expressions: 


g“=g*+%pt,  p'  = 


i 

y 

t 

L2(Q.) 

=pc(g*-j'(gt))- 


■g*» 


(3.8) 


where  PG  is  the  operator  of  projection  onto  the  admissible  convex  set  G.  For  a  rectangular  domain  Q  , 

the  projection  operator  periodically  cuts  those  values  going  beyond  Q  .  The  step  of  the  method  is 
chosen  in  accordance  with  the  rule 


a, 


■  mini 


"4/3 
+  £ 


g 


,  where  s  >  0  js  a  parameter  of  the  method. 


is 


If,  for  some  k,  gk =0,  then  the  necessary  condition  of  minimum,  gA  =Pc(gA  —  Jr(gA )) , 

achieved,  and  the  iterations  should  be  stopped.  The  following  statement  governs  the  convergence  of 
the  method. 

Theorem  3.4.  Let  the  conditions  of  theorem  3.3  hold  and,  additionally,  G  is  a  closed  convex  set 

k 

and  the  consequence  g  from  (3.3)  is  infinite.  Then 


pG(g‘  -  J'(g*))- 


L2(Q.) 


— »  0,  k  ->  oo . 


(3.9) 


Comment.  Condition  (3.10)  means  that  the  consequence  g'  from  (3.3)  converges  to  the 
asymptotic  version  of  the  necessary  condition  of  minimum  mentioned  above. 

For  practical  calculations,  we  also  used  a  simplified  gradient  projection  method,  in 
accordance  with  the  following  iteration  law: 

gk+l  =Pn(g*  ~ak3'(gk)).  (3  10) 

The  sequence  { ak }  of  steps  in  (3.10)  was  determined  by  the  splitting  method  on  the  basis  of  the 

condition  of  a  strict  decrease  of  the  functional  at  each  step:  «/( g^+l)<  J(  g*)- 

The  third  gradient-like  method  elaborated  is  a  conditional  gradient  method.  For  the  case  of 
the  rectangular  domain  Q  =  [0,2zr]  x  [0,2/r] ,  the  corresponding  iterative  sequence  is  constructed 
according  to  the  formula 

g£+l  (x)  =  g  k(*)  +  «(g  k(*)  ~  g/tW) »  (3  n) 

where  g^(x)  =  (gk  (x),gki  (x))  is  an  additional  vector  taken  from  the  following  cut-off  rule: 


8k  m  ^x) 


=  0,  J’m( g)>0, 

W,^w(g)<0, 


m  =  1,2  , 


J'm  (g)  being  the  m- th  coordinate  of  the  gradient  J'(g)  =  ( J\  (g),/'2  (g))  • 


2.3.5.  Projection  finite-element  approximations 
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In  the  scope  of  the  Project,  special  finite-dimensional  approximations  have  been  elaborated. 
During  the  investigation,  the  following  three  sufficient  obstacles  in  implementing  traditional  finite- 
difference  approximations  have  been  successfully  overcome: 

•  Superposition  with  an  argument  transform  in  (3.3)  acts  as  a  distribution  and  requires  an 
adequate  approximation  to  be  implemented; 

•  The  value  of  the  transformed  argument  g(x)  in  integral  (3.3)  may  coincide  with  none  of  the 
mesh  nodes; 

•  Even  for  an  infinitely  smooth  but  irreversible  argument  transform  (such  as  g  (x )  =  x  () 

that  is,  focusing  onto  a  point),  the  corresponding  solution  of  (3.1),  u(x,t ) ,  need  not  to  be  a 
smooth  enough  to  implement  classical  finite  differences. 

The  idea  of  the  approach  we  implemented  is  to  utilize  projection  finite-element  scheme,  which 
keeps  the  structure  of  distributed  kind  of  generalized  superposition  in  (3.1)  and  (3.3).  To  describe  it, 

we  use  the  space  S\  =  Lin{(pij(x)\  of  piece-wise  linear  functions  (pij(x)  =  (pi(xi)-(pj(x 2) 

constructed  on  the  basis  of  first-order  B-splines,  as  above.  The  approximate  solution  of  (3.1)  is  sought 
as  an  element  w  =  (w°,...,wM )  e  S  of  the  space  S  =  Sj  x5j  xSj  X...X.S)  ,  satisfying  the  following 

M+ 1  times 

finite-element  projection  analog  of  (3.1),  (3.3): 

(r~l(wm  -  pj  +  (wm ,  <pj  +  d(v j_wm ,  V_l^}  =  (F(wm ),  (p(g)j,  m  =  1 M  ,  V>  e  5j , 

where  brackets  mean  the  inner  product  in  the  Lebesgue  space  L  (Q) .  There  is  another  matrix  fonn  of 

the  previous  formula,  being  more  convenient  for  computations: 

{(\  +  T)B  +  rDh)wm  =  T$>(wm)  +  Bwm  \  m  =  .  (3.12) 

To  solve  nonlinear  matrix  equation.  (3.12),  we  applied  iterations: 

((l  +  r)5  +  rM)w(s+1)  =  r<D(ww)  +  5w"Hl,  5  =  0,1,...,  w(0)  =wm~l.  (3.13) 

It  is  important  to  underline  that  calculating  the  matrix  <b(/)  with  elements 


,l(f)  =  J  W(x))  n,l  (gW) ,  feS  1 » 

h\h2  J 


Q 

N\,N2 


leads  to  sums  ^  F(f(x™,x2  ))(Pkj(&(xT ’x2  ))  »  which  are  very  expensive  for  direct 

m,n=\ 


computations.  Basing  on  the  finiteness  of  spline  support,  substantial  economy  of  the  operations 
number  can  be  achieved.  Indeed,  since 

SUPP <Pk,l  =(xi~l,xi+l)x(xl2~l,xl2+l)  =  ^kj,  k  =  l  =  l,...,N2  -1, 

it  is  sufficient  to  sum  over  those  indexes  (m,n)  for  which  g(x™ ,x2)  eClkJ  and  to  keep  the 

information  of  prototypes  g\<Pk  /(D))  in  the  form  of  the  following  structure  (see  fig.  3.1): 

struct  fiGArray{ double  value;  struct  fiGArray*  next;  int  m,n;}, 
where  value  means  the  value  of  (p^  / (g( xf! , x2  )) ;  next  is  the  pointer  to  the  next  element;  m,  n  are 
the  node  numbers  before  the  transfonn. 


39 


Fig.  3.1.  Structure  for  keeping  the  prototype  information. 


In  the  scope  of  the  Project,  we  have  presented  a  theoretical  basis  for  the  projection  method  used. 
Theorem  3.5.  Let  the  conditions  of  Theorem  3.3  hold.  Then: 

1.  For  every  sufficiently  small  temporal  step  r>0,  Eq.  (3.12)  have  a  unique  solution.  The  iterations 

w(v)  given  by  (3.13)  converge  to  the  exact  solution  wm  of  Eq.  (3.12)  in  L2(Q),  according  to  a 

1/2 

geometric  progression  with  the  denominator  value  being  proportional  to  r  ; 

2.  The  estimate  of  convergence,  ||  u  —  w  Wji^Q  )  —  c(r14  a+\h  |1/_  2a),  is  valid  with  the  constant  C 


being  independent  of  r,  /?, ,  h2 . 

By  analogy  with  Eq.  (3.12),  a  projection-difference  scheme  has  been  suggested  and  elaborated 
for  numerically  solving  the  conjugate  problem  (3.6).  In  the  operator  form,  the  projection-difference 

scheme  for  searching  e  5)  can  be  expressed  as 


((1  +  t)B  +  rDA)c"‘-]  =  fm 


r = b(2  TP(tm  x  /  -  «r ) + Sm  ) + V 


(3.13) 


\T/  m  I  1  \p m  (  pt  m\ 

where  T  is  used  for  approximating  the  inner  product  \  '  v  j .  To  implement  the 

scheme,  it  is  essential  to  use  the  information  on  the  mutual  correspondence  between  the  prototype  of 
points  and  spline  supports,  which  is  built  during  solving  the  direct  problem  and  is  stored  in  the 
fiGArray  structure  described  above. 


2.3.6.  Testing  the  simplified  gradient  projection  method 

Algorithm  (3.10)  was  first  tested  with  spatially  ID  target  functions  for  the  problem  of  phase 
distortion  suppression,  described  by  Eq.  (3.1)  subject  to  In  -periodic  boundary  conditions.  Such  a 
distortion  was  introduced  to  a  spatially  inhomogeneous  (with  respect  to  argument  x\ )  input  phase 
profde: 

uq{xi,X2)  =  2  +  0.3(exp(-15(x1-7r)2)+exp(-15(x1-7r/3)2)+exp(-15(x1-57r/3)2  J. 

The  suppression  was  ensured  by  the  feedback  loop  of  a  nonlinear  interferometer  with  the  nonlinearity 
F(u )  =  A(1  + ;/  cos  a ) ,  where  A  is  a  coefficient  proportional  to  the  intensity  of  the  input  light  wave; 
and  y  is  the  contrast  of  the  interference  pattern.  In  the  computation  performed,  the  parameters  of  the 
nonlinearity  and  equation  (3.1)  were  set  to  K  -  3.6 ,  y  =  1 ,  D  =  0.07  .  It  was  necessary  to  find  the 
argument  transformation,  g(x) ,  which  ensures  the  minimum  deviation  of  the  phase  profile  u  =  u(x,  t ) 
from  its  initial  (at  t  -  0 )  distribution  liq (x)  within  the  whole  time  interval  t  e  [0,  T]  for  T  -  3  .  Such  a 
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statement  is  formulated  in  the  form  of  the  minimization  problem  (3.4)-(3.5)  with  the  stationary  target 
profile  u\(x,t)  =  uq(x)  .  The  weight  function  contained  in  (3.4)  had  the  following  form: 


p(x,t)  =  exp 


1  -T 


!( T 2  ~(t-T)2\ 


n 


which  takes  into  account  that  the  contribution  of  the  final  portion  of  the  trajectory  to  the  functional 
value  exceeds  that  of  the  initial  portion.  As  an  initial  approximation  for  the  control  law  in  the  gradient 

method  (15),  the  identical  transform,  g°(x)  =  x,  was  taken.  By  the  final  point  of  time,  the  solution 


became  spatially  homogeneous,  with  the  functional  J( g°  )  =  0.79  . 

As  a  result  of  two  steps  of  the  gradient  method  (3.14),  the  value  of  the  functional  decreased 

down  to  J{ g2  )  =  0.30,  that  is,  by  the  factor  of  more  than  2.5.  Figs.  3.2a-b  show  the  profiles  of  the 
target  function  and  the  optimized  solution  obtained.  Cross-sections  of  these  profiles  are  shown  in  Fig. 
3.3a.  It  can  be  seen  that  the  optimization  procedure  has  succeed  in  substantially  reflecting  the  main 
features  of  the  target  function.  One  can  notice  certain  symmetry  peculiar  to  the  optimization  problem: 
if  the  initial  approximation  of  the  transform  does  not  differ  from  the  identical  with  respect  to  variable 
X2  and  if  the  target  profile  is  independent  of  this  variable,  then  the  solution  of  the  conjugate  problem 
(and,  consequently,  the  functional  gradient  as  well)  is  also  independent  of  X2 .  Therefore,  the  gradient 
procedure  (3.14)  generates  a  sequence  of  argument  transforms  for  which  the  first  component,  g\  (x) ,  is 
responsible  for  spatial  nonhomogeneity  of  the  target  profile,  while  the  second  component  does  not 
change:  g2(x)  =  x2  •  Figs.  3.2c-d  show  the  profiles  of  the  deviations  of  the  corresponding  components, 

g]  (x)  and  g2(x) ,  of  the  second  iteration  of  the  argument  transform,  g  (x)  =  (g1  (x),  g2(x)) ,  from  the 

identical  transform,  g°(x)  =  x .  The  cross-section  of  gi(x)  is  shown  in  Fig.  3.3b. 


a.  Phase  profile  m(x,J).  c.  Profile  gi(x)-xi. 


b.  Target  profile  u\(x,T).  d.  Profile  g2(x)-X2- 


Fig.  3.2.  Results  of  optimization  for  the  problem  of  phase  distortion  suppression. 
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Fig.  3.3.  Cross-sections  of  x2  =  const:  u(x,T)  and  ui(x,T )  from  Fig.2a-b  (a); 
gj(x)  component  of  the  optimized  transform  (b). 


In  the  second  example,  the  control  target  was  switching  system  (3.1)  from  the  initial  spatially 
homogeneous  state  u(x,t  =  0)  =  uq(x)  =  2  with  a  constant  (for  each  point  x)  velocity  to  the  final  (at 

T  -2)  state  defined  by  a  phase  profile  having  four  localized  peaks  (see  Fig. 3. 4).  In  the  corresponding 
functional  (3.4),  the  target  state  was  defined  by  the  following  formula: 

U](x],x2)  = 

=  2  +  t-  (exp(-  9(x\  - 2tt/3)2  -  9(x2  -  2tt/3)2)+  exp(-  9(x:  -  4tt/3)2  -  9(x2  -  2tt/3)2 ))+ 

+  t  •  (exp(-  9(x!  -  2tt/3)2  -  9(x2  -  4tt/3)2  )+  exp(-  9(xj  -  4tt/3)2  -  9(x2  -  4tc/3)2  )) 

In  the  case  of  the  identical  transformation  g°(x)  =  x,  the  solution  of  problem  (3.1)  remained  spatially 
homogeneous  within  the  whole  time  interval  t  e  [0,  T] ,  which  led  to  the  functional  value  of 

J(g° )  =  1 .63  (for  the  following  parameters:  K  -  4.3 ,  y  =  1 ,  D  -  0.07 ).  As  a  result  of  five  steps  of  the 

gradient  method  (3.14),  the  functional  decreased  down  to  /(g5  )  =  0.37  ,  that  is,  more  than  4  times. 
The  optimized  solution  profile  u(x,T)  well  corresponded  with  the  target  profile  u\  (x,  T)  (see 

Fig.3.4a-b).  The  deviations  of  the  obtained  components,  g|  (x)  and  g2(x),  of  the  transform 

g5(x)  =  (g1  (x),g2(x))  are  shown  in  Fig.3.4c-d.  It  can  be  seen  that  the  x1  -cross-section  of  the  g1  (x) 

profile  and  the  x2  -cross-section  of  the  g2(x)  profile  look  very  similar  to  those  plotted  in  Fig. 3. 3, 
especially  at  their  localization  peaks. 
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a.  Phase  profile  «(x, 7). 


c.  Profile  gi(x)-X|. 


b.  Target  profile  u\(x,T). 


d.  Profile  g2(\)-x2. 


Fig.  3.4.  Optimization  results  for  a  2-D  localized  target  function. 

Basing  on  the  above  calculations,  the  following  conclusions  can  be  made  on  the  features  of 
formation  of  localized  solutions  by  means  of  the  spatial  argument  transform  control: 

3.  Localized  structures  can  be  obtained  by  means  of  argument  transforms  that  are  not  invertible  in 
the  vicinity  of  each  localization  peak.  A  similar  conclusion  has  been  earlier  made  for  the 
spatially  homogeneous  case. 

4.  If  the  initial  and  the  target  state  of  the  system  are  independent  of  one  of  the  variables,  then  the 
corresponding  component  of  the  optimized  transform  can  be  considered  identical  with  this 
variable. 

5.  Simultaneous  localization  with  respect  to  two  variables  leads  to  an  optimized  transform,  each 
component  of  which  is  mainly  responsible  for  the  localization  along  the  corresponding 
variable. 

2.3.7.  Testing  the  conditional  gradient  method  in  contrast  with  the  simplified  method 

We  have  tested  and  compared  the  gradient  projection  and  the  conditional  gradient  method  for 
the  problem  of  constructing  a  phase  profile  being  maximally  close  to  the  desired  target  profile  at  the 
final  time  point  T  .  For  this  model,  it  is  convenient  to  use  as  a  weight  function  p(t  )  some  function  that 
is  nonzero  only  within  a  portion  of  the  segment  and  equals  unity  at  the  final  time  point.  Steeper 
functions,  such  as  super-Gaussian,  are  more  advantageous  for  this  purpose.  When  modeling  the 
optimization  problem  with  the  following  parameters 

ui(x\,X2,t)  =  2.02+0.25t(x^2+X22)3  exp(-2  •  (x^2+X22)4)  , 

D  =  0.02,  K  =  3.6,  /  =  I 

uq{x\,X2)  =  2.02  +  0.5(xi2  +  X22)3exp(-2  •  (x^2  +  X22)4) 
g°(x)  =  x,  T  =  2, 

by  means  of  the  gradient  projection  method,  the  target  functional  finally  decreased  by  90.8%.  When 
using  the  conditional  gradient  method,  the  target  functional  decreased  by  70.9%. 
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Fig.3.5.  Optimization  results  for  an  annular  target  profile:  target  function  (a);  solution  and 
transform  components  obtained  by  the  gradient  projection  method  (bl-b3);  solution  and 
transform  components  obtained  by  the  conditional  gradient  method  (cl-c3). 

We  also  studied  possibilities  to  model  complex  structures  containing  both  rectangular  and 
circular  regions.  In  this  case,  the  conditional  gradient  method  shows  better  results.  For  the  example 
described  below,  this  method  yielded  a  65.2%  target  functional  decrease,  while  the  gradient  projection 
method  yielded  61.5%. 
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Fig.3.6.  Optimization  results  for  a  target  profile  consisting  of  a  rectangular  and  a  circular 
region:  target  function  (a);  solution  and  transform  components  obtained  by  the  gradient 
projection  method  (bl-b3);  solution  and  transform  components  obtained  by  the  conditional 
gradient  method  (cl-c3). 

The  parameters  used  for  the  simulations,  the  results  of  which  are  shown  in  Fig.3.6,  were  the 
following: 

D  =  0.01,  K  =  3.6,  y  =  0.8, 

uq(xi,X2)  =  2  +  exp(-50(xj2  +(0.1x2  -^/3)2)(xi2  +(*2  -^/3)2)4)  +  exp(-18xi6  -5(x2  +  zr/5)6) 
g°(x)  =  x,  T  =  2, 

u\(xi,X2,t)  =  2  +  0.5t(exp(-50(x^2  +(0.1x2  -^/3)2)(xi2  +  (*2  -^/3)2)4)  +  exp(-18x^6  -5(x2  +  zr/5)6))  . 

Although  the  solution  obtained  by  the  conditional  gradient  method  better  reproduces  the  shape  of  the 
target  profile,  this  solution  is  well  perturbed  (like  in  the  previous  example). 

2.3.8.  Testing  special  gradient  projection  method  in  contrast  with  the  simplified  method 

Example  1.  In  the  numerical  simulation  series  discussed  below,  the  target  function  was 
defined  by  the  following  formula: 

ux{xvx2,t)  =  2  +  exp(-((xj  - 7t)2  +(x2  -^•)2))cos((x1  - n) - (x2  -rc){{xx  -re)1  +(x,  -7t)2)2). 

A  plot  of  the  target  phase  distribution  for  the  moment  t  =  T  is  shown  in  Fig.  3.7. 
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Fig.  3.7.  Target  function  ut  (xvx2,t  =  T ) 

The  parameters  used  for  the  simulations,  the  results  of  which  are  shown  in  Fig. 3. 8,  were  the 
following: 

ul(xl,x2,t)  =  2  +  exp(-((xj  -n')1  +  (x,  -^■)2))cos((x1  - n) - (x2  -7i)((xx  -n')2  +  (x2  -7i)2)2), 
w0(x1,x2)  =  2,  D  =  0.01,  K  =  3.2,  y  =  0.7,  g°(x)  =  x,T  =  4. 
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Fig.  3.8.  Optimization  results:  simplified  projection  gradient  method  (al-a3)  and  special 
projection  gradient  method  (bl-b3). 

In  this  series  of  numerical  simulations,  both  methods  showed  similar  results  at  different  values 
of  the  step  splitting  factor  of  the  gradient  method  (the  system  parameters  being  the  same).  After  four 
steps  of  the  gradient  method,  the  classical  gradient  projection  method  yielded  a  71,48%  decrease  of 
the,  while  the  modified  method  ensured  a  71,65%  decrease.  As  it  can  be  seen  from  Fig.  3.8,  the  results 
are  hardly  distinguishable  visually. 

Example  2.  In  the  next  series  of  numerical  simulations,  the  target  function  was  defined  as 
ul(xl,x2,t)  =  2.2+0.22 sin((xj  -n)  +  1.4cos(3(x2  —  x2 ))) .  A  plot  of  the  target  phase  profile  is  shown  in 

Fig.3.9. 
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Fig.3.9.  Target  function  ut  (xvx2,t  =  T ) 


The  parameters  used  for  the  simulations,  the  results  of  which  are  shown  in  Fig. 3. 10,  were  the 
following: 

D  =  0.05,  K  =  4,  y  =  1,  w0(xl5x2)  =  2.2+0. 22 sin((x1  -x)  + 1.4 cos(3(x2  —  jcx ))),  g°(x)  =  x,  T  =  4 
u1(xl,x2,t)  =  2.2+0.22  sin^Xj  -;r)  +  1.4cos(3(x2  -x,))). 


Fig.  3.10.  Optimization  results:  simplified  projection  gradient  method  (al-a3)  and  special 
projection  gradient  method  (bl-b3). 

After  9  iterations,  the  simplified  method  yielded  a  52,98%  decrease  of  the  target  functional  with 
respect  to  its  value  before  the  optimization.  At  the  same  time,  the  special  method  was  capable  of 
decreasing  the  functional  by  only  by  40,79%,  after  3  iterations.  Besides,  in  this  case  the  special 
method  appeared  to  be  less  flexible:  with  its  use,  the  optimization  results  remained  almost  the  same 
regardless  of  the  step  splitting,  while  with  the  simplified  method  allowed  us  to  improve  the 
optimization  results  by  step  splitting. 

Example  3.  In  this  series  of  numerical  simulations,  the  target  function  was  defined  as 
ul(xl,x2,t)  =  1.8  +  0. 3exp(sin((x1  -k){x2  -^•)(5(x1  -7i)~ 4(x2  -^r))  +  (xl  -n')1  +(x2  —  zr)2 )  —  1). 

A  plot  of  the  target  phase  profile  is  shown  in  Fig.3.11: 
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Fig.3.11.  Target  function  ul(xvx2,t  =  T ) 

The  parameters  used  for  the  simulations,  the  results  of  which  are  shown  in  Fig. 3. 12,  were  the 
following: 

D  =  0.02,  K  =  3.2,  y  =  1,  u0(xux2)  =  2,  g°(x)  =  x,  T  =  4, 

=  1.8  +  0. 3exp(sin((xj  -^)(jc2  -^)(5(Xj  -^)-4(jc2  -^))  +  (jCj  -;r)2  +(x2  -;r)2)-l). 


Fig.  3.12.  Optimization  results:  simplified  projection  gradient  method  (al-a3)  and  special 

projection  gradient  method  (bl-b3). 

In  this  series  of  numerical  simulations,  the  special  method  shoed  better  perfonnance,  yielding  a 
41,84%  decrease  of  the  target  functional  after  only  2  iterations.  The  same  2  iterations  performed  for 
the  simplified  method  yielded  only  a  21,27%  decrease  of  the  functional.  Visually  comparing  the 
results,  one  may  notice  that  the  special  method  provides  better  resolution  of  fine  details  (see 
Fig.3.12bl),  while  the  simplified  method  blurred  these  (see  Fig.3.12al).  . 

The  results  of  the  above  simulations  demonstrate  that  no  universal  gradient-type  method 
obviously  overscoring  the  other  methods  does  exist.  The  special  gradient  projection  method  is 
theoretically  convergent  regardless  of  the  initial  iteration  and  is  more  formalized.  However,  it  has 
fewer  fitting  parameters  and  therefore  is  less  flexible.  On  the  contrary,  no  convergence  is  theoretically 
guaranteed  for  the  simplified  method.  However,  the  wide  possibilities  of  step  splitting  allow  one  to 
achieve  good  results  with  this  method. 

1.0.0.  Static  and  iterative  Tikhonov  regularization  methods. 
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In  the  studies  of  the  conditional  gradient  method,  the  shape  of  the  solutions  was  perturbed, 
especially  within  constant  areas.  To  obtain  more  smooth  solutions,  we  propose  to  employ 
regularization.  To  investigate  its  capabilities,  we  used  a  regularized  version  of  the  conditional  gradient 
method.  The  functional  to  minimize  is  the  Tikhonov  functional 


Tk(g)  =  J(g)  +  ak  ||g|f2  ->mf , 

^  G 


(3.14) 


while  for  its  gradient  the  following  formula  is  valid:  7)  ’(g)  =  ./'(g)  +  2a k g ,  where  the  target  functional 
gradient  is  calculated  according  to  (3.7)  with  the  nonlinearity  F(u)  =  K{\  +  y  cosu) ,  which 
corresponded  to  the  nonlinear  interferometer  model,  while  the  computation  time  interval  was  T  —  1.5 
(in  the  units  of  the  response  time  of  the  nonlinear  medium). 

First,  we  employed  static  regularization  with  the  regularization  parameter  ak  being  constant, 

regardless  of  the  iteration  number.  These  studies  showed  that  the  solution  obtained  with  the 
regularization  method  has  a  more  smooth  profile.  However,  such  regularization  reduced  the 
“flexibility”  of  the  gradient  method,  that  is,  increased  the  number  of  the  steps  required  and  decreased 
the  relative  reduction  of  the  functional  (relatively  to  the  case  of  no  regularization). 
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Fig.  3.13.  Solution  and  transform  components  obtained  by  the  conditional  gradient  method 
without  regularization  (al-a3)  and  with  regularization  (bl-b3). 

The  parameters  used  for  the  simulations,  the  results  of  which  are  shown  in  Fig. 3. 13,  were  the 
following: 

D  =  0.02,  K  =  4.2,  /  =  0.9,  uq(xi,x2)  =  2.2  , 

u1(xl,x2,t)  =  2.2  +  0.8(exp(-15(x,4  +  x24  +  4xfix.fi)2)  +  exp(-5(xj4  +  x24  +  4xfixfi)2) . 

In  the  case  of  no  regularization,  14  steps  of  the  gradient  method  were  made,  which  led  to  27% 
reduction  of  the  target  functional.  When  we  “turn  on”  the  regularization,  the  functional  value 
decreased  by  23.75%  after  9  steps.  Such  a  regularity  can  be  explained  by  the  fact  that  the  functional 
contains  an  additional  summand  whose  contribution  to  the  whole  functional  increases  as  the  solution 
approaches  the  target  function.  This  reduces  the  possibilities  of  perfectly  matching  the  solution  with 
the  target  profile. 

Some  numerical  studies  were  performed  in  order  to  reveal  the  influence  of  this  summand 
(namely  -  the  regularization  parameter)  on  the  matching  quality.  The  results  of  these  studies  showed 
that  increasing  the  regularization  parameter  decreased  both  the  number  of  the  gradient  method  steps 
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and  the  matching  quality.  The  latter  fact  shows  up  in  the  relative  decrease  of  the  functional  value,  as  it 
is  illustrated  in  the  table  below. 


Functional  decreasing,  %  Regularization  parameter 

Number  of  steps 

32,16 

0,00005 

7 

30,94 

0,0001 

6 

26,07 

0,0005 

4 

Tab.3.1.  Influence  of  the  static  regularization  parameter  on  the  optimization  capabilities. 

The  results  contained  in  this  table  correspond  to  the  following  parameters  of  the  model: 

D  =  0.007,  K  =  4,  y  =  1, 

u0(x  j,x2)  =  2  + 

+  exp(-9(xj2  +(0.1x,2  +7i/7)2(xl2  +(x22  +n!  7)2(x,~  +(x,2  +  ;z/7)2)(x12  +(x22  +  7ill)2)))) 
ul(x1,x1,t)  =  2  + 

+  exp(-9(X[2  +(0.1x22  +^/7)2(Xj2  +(x22  +  zz/7)2  (Xj2  +(x22  +  7i/7)2)(x2  +(x22  +  ;z/7)2)))) 

To  solve  the  optimization  problem  more  efficiently,  iterative  regularization  was  employed. 
As  the  iterative  sequence,  we  chose  the  following  one:  ak  =  a0(k  + 1)_/  ,  where  k  is  the  number  of  the 
current  step  of  the  gradient  method,  while  a0  and  y  are  parameters  of  the  numerical  experiment. 

Thus,  the  influence  of  the  additional  summand  decreases  as  the  step  number  (iteration  number) 
increases.  The  below  example,  which  corresponds  to  the  following  values  of  the  parameters 

D  =  0.02,  K  =  3.6,  y  =  l,  uq{x\,X2)  =  2.02  +  (x^2  +  X22)3exp(-2  •  (x^2  +  x^)^) 

ul(xl,x2,t)  =  2.02+0.5?(x1"+x22)3  exp(-2  •  (X[“+x22)4) , 
allows  one  to  compare  the  results  of  numerical  experiments  involving  no  regularization,  static 
regularization,  and  iterative  regularization. 
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Fig.3.14.  Optimization  results  for  an  annular  target  profile:  target  function  (a);  solution  and 
transform  components  obtained  by  the  gradient  projection  method  without  regularization 
(bl-b3);  using  static  regularization  (cl-c3);  using  iterative  regularization  (dl-d3). 

In  the  case  of  static  regularization,  the  functional  decreased  by  26.4%;  in  the  case  of  iterative 
regularization  -  by  50%;  in  the  case  of  no  regularization  -  by  70.9%. 

Thus,  in  the  framework  of  this  direction  of  investigations,  we  studied  the  regularization  method 
aimed  at  obtaining  a  more  smooth  solution.  For  this  purpose,  methods  of  both  static  and  iterative 
regularization  were  employed.  Our  studies  showed  that  although  static  regularization  yielded  a 
smoother  solution,  the  presence  of  the  addition  summand  in  the  target  functional  somewhat  reduced 
the  optimization  capabilities.  In  our  opinion,  the  best  alternative  to  the  static  regularization  is  iterative 
regularization  that  both  allows  one  to  reasonably  control  the  contribution  of  the  additional  summand 
and  yields  a  smoothed  solution  that  is  more  close  to  the  target  function  than  the  solution  obtained  with 
static  regularization. 

These  studies  showed  that  when  we  solved  the  minimization  problem  using  the  gradient 
projection  method,  the  structure  of  the  argument  transform  being  the  optimization  result  was  close  to 
the  structure  of  the  target  function.  When  the  conditional  gradient  method  was  used,  the  resultant 
argument  transform  was  not  too  similar  to  the  target  profile  structure.  This  was  due  to  the  specific 
character  of  forming  a  new  approximation  in  the  conditional  gradient  method.  Thus,  although  the 
conditional  gradient  method  yields  more  diverse  possibilities  to  perform  the  optimization,  the  target 
profile  structure  is  reproduced  more  precisely  when  using  the  gradient  projection  method.  To  join  the 
potentials  of  the  two  methods,  we  proposed  a  new  problem  statement,  in  which  some  approximation  of 

the  argument  transform,  g ,  chosen  on  the  basis  of  a  priori  information,  participated  in  the  stabilizing 
term  ||  g  -  g  ||2  as  an  additional  parameter.  Thus,  we  suggest  replacing  functional  (3.14)  with 

Tk(g)  =  J(g)  +  ak  ||  g  -g  1 1 2  > inf ,  (3.15) 

( j 

where  g  is  an  approximation  of  the  argument  transform,  chosen  on  the  basis  of  a  priori  information. 
Our  primary  goal  was  to  figure  out  how  this  parameter  affects  the  optimization  quality. 

The  influence  of  the  parameter  g  on  the  optimization  quality  was  studied  numerically.  The 
corresponding  numerical  simulations  involved  a  few  stages.  First,  the  problem  was  solved  using  the 
gradient  projection  method  with  no  regularization.  Then,  the  structure  of  the  argument  transform 
obtained  was  analyzed  in  order  to  build  the  approximation  g  being  the  parameter  of  the  stabilizing 
term.  Afterwards,  optimal  control  problem  (3.15)  was  solved  using  the  conditional  gradient  method  for 
the  case  of  iterative  regularization. 
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Example  1.  The  parameters  used  for  the  simulations,  the  results  of  which  are  shown  in 
Fig. 3. 15,  were  the  following:  g°(x)  =  x,r  =  2. 

D  =  0.02,  K  =  3.6,  7  =  1, 

w0(x j,x2)  =  2.02  +  ((X[  -;r)2  +  (x2  -;r)2)3exp(-2  •  ((xj  -  n'f  +  (x2  -;r)2)4) 
ufxvx2,t)  =  2.02+0.5t((x1  -;r)2  +  (x2  -7r)2)3exp(-2-  ((Xj  -  n'f  +(x2  -;r)2)4) 
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Fig.3.15.  Optimization  results  for  an  annular  target  profile,  obtained  by  the  conditional 
gradient  method:  solution  and  transfonn  components  obtained  with  the  additional  parameter 
contained  in  the  stabilization  term  (al-a3)  and  without  the  additional  parameter  (bl-b3). 

When  we  used  no  additional  parameter  in  the  stabilizing  term,  the  target  functional  decreased  by 
50.1%  after  7  steps  of  the  conditional  gradient  procedure.  When  the  additional  parameter  was 
“switched  on”,  the  functional  decreased  by  74.65%  after  only  3  steps.  The  following  approximation 
was  taken  as  the  additional  parameter: 

g(Xj)  =  Xj  +  0.2sin(x1)((x1  -nf  +  (x2  -;r)2)3exp(-2  •  ((Xj  -nf  +  (x2  -;r)2)4) 

g(x2)  =  x2  +  0.2sin(x2)((Xj  -nf  +(x2  -^')2)3exp(-2-((x1  -nf  +(x2  -;r)2)4) 

Example  2.  g°(x)  =  x,  T  =  1.5 ,  D  =  0.02,  if  =  4.2,  y  =  0.9,  w0(x1,x2)  =  2.2 
w1(x1,x2,t)  =  2.2  +  0.8(exp(-15((x1  -nf  +  (x2  -nf  +  4(xj  -;r)2(x2  -;r)2)2)  + 

+  exp(-5((xj  -  nf  +  (x2  -  nf  +  4(x,  -nf{x2  -nff). 
g(Xj)  =  Xj  +0.2sin(x1)(exp(-15((x1  -  nf  +(x2  -  nf  +  4(xj  -nf(x2  -  nff)  + 

+  exp(-5((x1  -  nf  +(x2  -  nf  +  4(xj  -;r)2(x2  -nff) 
g(x2)  =  x2  +0.2sin(x2)(exp(-15((x1  -  nf  +  (x2  -  nf  +  4(Xj  -;r)2(x2  -;r)2)2)  + 

+  exp(-5((xj  -  nf  +  (x2  -  nf  +  4(X[  -;r)2(x2  -nff) 
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Fig.3.16.  Solution  and  transform  components  obtained  by  the  conditional  gradient  method 
with  regularization  using  the  additional  parameter  in  the  stabilization  term  (al-a3)  and 
using  no  additional  parameter  (bl-b3). 

Without  the  additional  parameter,  the  target  functional  decreased  by  23.75%  after  9  steps  of  the 
conditional  gradient  procedure.  When  the  additional  parameter  was  “switched  on”,  the  functional 
decreased  by  35.8  %  after  8  steps. 

In  both  numerical  experiments,  a  quite  simple  approximation  of  the  argument  transform  was 
used:  if  ul(x,t)  =  C  +  f(t)u(\) ,  where  C  is  a  constant  and  fit)  is  a  time-dependent  function,  then 
g(x)  =  X  +  Asin(x)w(x) ,  where  the  amplitude  A  is  chosen  using  the  information  on  the  transfonn, 
obtained  with  the  gradient  projection  method. 

Thus,  for  the  case  of  localized  optical  patterns  with  a  fine  structure  of  the  solution  profile,  the 
proposed  method  allows  one  to  obtain  better  result  than  that  obtained  with  the  conditional  gradient 
method,  both  with  regularization  and  without  it.  Moreover,  the  proposed  method  requires  fewer 
iterations.  A  quite  simple  approximation  can  be  taken  as  the  additional  parameter  g .  So  far,  we 
haven’t  obtained  comparably  good  results  for  spatially  distributed  structures  and  localized  structures 
having  large  spatially  uniform  areas.  This  is  because  a  more  careful  selection  of  the  approximation  is 
required  in  this  case. 
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3.  Results  and  Conclusions 


3.1.  In  the  scope  of  Task  1,  the  original  projection  method  of  integral  Fourier  transformation 
based  on  expanding  the  solution  into  a  series  of  a  full  system  of  Fourier  transform  eigenfunctions 
(Hermite  functions)  shows  very  promising  results  in  image  filtering,  texture  discrimination,  and  noisy 
image  edge  detection.  The  joint  use  of  this  method  with  Tikhonov  regularization  algorithms  aloows  us 
to  construct  very  stable  methods  for  image  processing  and  analysis. 

3.2.  In  the  scope  of  Task  2,  a  profound  mathematical  study  has  been  accomplished  for  the 
problem  of  optimal  discrete  Fourier  filtering.  It  covers  novel  results  on  the  existence  and  uniqueness 
theorems  for  both  the  direct  and  the  conjugate  initial-boundary  value  problems  for  the  functional 
differential  parabolic  equation.  The  solvability  of  the  filter  optimization  problem  has  been  proved  for  a 
Hilbert  brick  and  for  a  ball  of  the  radius  R  >  0  in  the  1 2  space  as  admissible  filter  sets.  The  terminal 
functional  gradient  formula  has  been  obtained  along  with  the  estimation  of  the  remainder  term.  The 
projection  finite-element  scheme  has  been  developed  both  for  direct  and  conjugate  problems  in  order 
to  figure  out  the  perfonnance  of  well-known  Fourier  filters  in  optical  systems  aimed  at  suppressing 
phase  distortions  of  the  input  light  wave.  Conditional  gradient  iterative  procedure  has  been  developed 
and  applied  for  searching  for  the  optimal  Fourier  filter  to  control  the  feedback  optical  system. 

3.3.  In  the  scope  of  Task  3,  we  have  suggested  and  developed  the  mathematical  statement  of 
the  optimal  control  problem  for  a  distributed  spatial  argument  transform.  It  is  applied  to  study 
controlling  optical  systems  with  nonlocal  transforms  of  a  light  wave  in  a  feedback  loop.  The  approach 
developed  utilizes  a  generalized  way  of  detennining  the  argument  transform.  As  an  advantage,  one  can 
use  a  wide  range  of  nonsmooth  and  irreversible  argument  transfonns,  and  consequently  achieve  better 
results  with  the  optimal  control.  Moreover,  this  approach  allows  us  to  develop  projection  finite- 
element  methods  for  approximating  both  the  direct  and  the  conjugate  problem  in  a  similar  way.  As  a 
result,  we  have  elaborated  and  theoretically  analyzed  computational  versions  of  the  projection  gradient 
and  the  conditional  gradient  method  for  target  functional  minimization. 

6.  Future  Work  Recommended 

The  efficiency  of  the  original  projection  method  of  integral  Fourier-Hermite  transformation  with 
Tikhonov  regularization  capability  continued  increasing  as  we  proceed  with  the  Project.  Nevertheless, 
even  the  speed  of  calculations  for  the  projection  method  was  optimized  during  the  Project,  we  see  a 
principal  possibility  to  construct  a  faster  version  of  the  algorithm.  The  idea  is  reducing  the  volume  of 
calculations  using  a  special  grid  for  expansion  coefficients  calculation.  The  grid  is  to  be  formed  using 
zeroes  of  Hennite  functions,  stored  as  an  additional  array.  Gibbs  (ringing)  effect  is  the  next  unpleasant 
property  of  the  oscillatory  sets  of  basic  functions.  Construction  of  Gibbs  effect  suppression  algorithms 
using  the  Tikhonov  regularization  method  and  other  postprocessing  methods  will  allow  us  to  enhance 
the  projection  method. 

The  preliminary  modeling  of  the  system  with  an  additional  integral  feedback  loop,  perfonned  in  the 
scope  of  the  Project,  showed  efficiency  of  this  system  for  the  task  of  phase  distortion  suppression.  The 
combined  use  of  the  integral  feedback,  Fourier  filtering  control,  and  argument  transfonn,  which  were 
development  in  the  framework  of  this  Project,  will  allow  one  to  considerably  improve  the  quality  of 
phase  distortion  suppression  and  beam  profile  shaping. 

Huge  potentials  of  improving  the  efficiency  of  the  operation  of  optical  feedback  systems  can  be 
revealed  by  studying  the  time-dependent,  continuously-corrected  control,  which  purposefully  correct 
the  dynamics  of  the  system  in  order  to  improve  the  control  quality.  Implementing  such  an  approach, 
which  is  widely  used  in  adaptive  optics,  with  the  use  of  time-dependent  Fourier  filters  and  nonlocal 
transfonns  of  spatial  arguments  will  allow  considerably  increasing  the  performance  of  the  system  and 
will  require  special  theoretical  and  numerical  studies. 


54 


5.  Personnel  Supported. 


Razgulin 

Alexander 

Vitalievich 

Project  director,  Researcher 

Krylov 

Andrey 

Serdjevich 

Researcher 

Denisov 

Alexander 

Mikhailovich 

Researcher 

Vasiliev 

Fedor 

Pavlovich 

Researcher 

Potapov 

Mikhail 

Mikhailovich 

Researcher 

Larichev 

Andrey 

Viktorovich 

Researcher 

Nikolaev 

Ilia 

Petrovich 

Researcher 

Chushkin 

Vladimir 

Alexandrovich 

Researcher 

Savvina 

Svetlana 

Sergeevna 

Researcher 

Kutovoi 

Andrey 

Vyacheslavovich 

Student 

Tsibanov 

Vladimir 

Nikolaevich 

Student 

5.  Technical  Publications. 

Papers: 

1.  A.V.  Razgulin  and  S.S.  Savvina.  Numerical  optimization  of  the  two-dimensional 
transfonnation  of  arguments  in  the  functional-differential  diffusion  equation  //  Computational 
Mathematics  and  Modeling.  2004.  Vol.  15.  No  4.  P.  333-343. 

2.  A.M.  Denisov,  A.S.  Krylov.,  and  V.N.  Tsibanov.  Edge  detection  method  by  Tikhonov 
regularization  //  Proceedings  of  14th  International  Conference  GraphiCon’2004.  P.  163-166. 

3.  A.V.  Razgulin  and  V.A.  Chushkin.  On  the  optimal  Fourier  filtration  for  a  class  of  models  of 
nonlinear  optical  systems  with  feedback  //  Computational  Mathematics  and  Mathematical 
Physics.  2004.  V.  44.  No.  9.  P.  1528-1538. 

4.  V.A.  Chushkin.  On  the  solvability  of  the  problem  of  Fourier  filter  optimization  for  one  class  of 
models  of  nonlinear  optical  systems  with  feedback  //  In  «Modeling  and  analysis  in  decision¬ 
making  problems)).  Moscow:  A.A.  Dorodnitsin  Computation  Center  of  RAS,  2004.  P.  123-135. 
(In  Russian). 

5.  V.A.  Chushkin.  On  the  attractor  of  the  equation  of  optical  Fourier-filtration  //  Vestnik  Mosk. 
Universiteta.,  ser.  15,  Computational  Mathematics  and  Cybernetics.  2005.  No  2.  P.  16-25  (In 
Russian). 

6.  A.V.  Razgulin.  On  parabolic  functional-differential  equations  with  controlled  transforms  of 
spatial  arguments  //  Russian  Academy  Doklady.  2005.  Vol.  403.  No.  4  (In  Russian). 

7.  A.V.  Razgulin.  The  control  problem  of  2D  transforms  of  spatial  arguments  for  parabolic 
functional-differential  equations  //  Differential  Equations  (In  Russian,  submitted). 

8.  A.V.  Razgulin.  Projection-difference  scheme  for  parabolic  functional  differential  equation  with 
2D  transform  of  arguments  //  Computational  Mathematics  and  Mathematical  Physics 
(submitted). 

9.  A.V.  Razgulin  On  a  gradient  projection  method  for  quasi-differentiable  functionals  with 
Holder-continuous  gradient  //  Moscow  University  Computational  Mathematics  and  Cybernetics 
(submitted). 


Abstracts: 

1 .  V.A.  Chushkin  On  optimization  of  a  discrete  Fourier  filter  //  Abstracts  of  XXVI  Conference  of 
Young  Scientists  of  Mechanics  and  Mathematics  Faculty  of  MSU,  Moscow,  2004.  P.  140  (in 
Russian). 

2.  V.A.  Chushkin  and  A.V.  Razgulin.  On  the  problem  of  Fourier  filter  optimization  in  a 
functional-differential  diffusion  equation  //  Modern  Methods  for  the  Theory  of  Boundary  Value 


55 


Problems.  Materials  of  Voronezh  Spring  School  «Pontryagin  Readings-XV».  Voronezh,  VSU, 
2004.  P.  235  (in  Russian). 

3.  V.A.  Chushkin.  Optimal  Control  Problem  by  Discrete  Fourier  Filter  //  Abstracts  of  4-th 
International  conference  on  operations  research.  Moscow,  September  21-24,  2004.  Moscow: 
MAKS  Press,  2004.  P.  54-58. 

4.  V.A.  Chushkin.  On  numerical  optimization  of  a  Fourier  filter  in  a  functional-differential 
equation  of  nonlinear  optics  //  Abstracts  of  XX  All-Russian  Workshop  « Analytical  methods 
and  optimization  of  processes  in  gas  and  liquid  mechanics  »,  Abrau-Durso,  September  4-7, 
2004.  Novosibirsk:  M.A.  Lavrentyev  Hydrodinamics  Institute  of  SB  RAS,  2004.  P.  79.  (In 
Russian). 

5.  V.A.  Chushkin,  A.V.  Larichev.,  I.P.  Nikolaev,  and  A.V.  Razgulin.  On  application  of  Fourier 
filtering  for  forming  a  given  phase  distribution  in  nonlinear  optical  systems  with  feedback  // 
Abstracts  of  the  III  International  Conference  «Basic  Problems  of  Optics  ».  St.  Petersburg. 
October  18-21,  2004.  St.  Petersburg:  ITMO,  2004.  P.  62-64.  (In  Russian) 


6.  Interactions/Transitions 

One  report  on  National  University  of  Singapore  -  MSU  forum  (December,  2003) 

A.  Razgulin.  Nonlinear  functional-differential  parabolic  equations  with  transfonned  spatial 
arguments:  stability,  control  and  finite-dimensional  approximations 

7.  Patent  Disclosures.  N/A 

8.  Honors.  N/A 


56 


